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ABSTRACT 


The  purpose  of  this  research  is  to  design,  simulate  and 
implement  a  robust  autopilot  system  for  the  vertical  mode  of 
operation  of  the  Archytas  prototype.  Archytas  is  an  Unmanned 
Air  Vehicle  that  is  designed  to  take  off  and  land  vertically, 
and  to  transition  to  horizontal  forward  flight.  A  feedback 
control  scheme  is  designed  for  both  the  single-input,  single¬ 
output  and  the  multi-input,  multi-output  subsystems  using 
optimal  control  techniques.  In  this  research,  the  linear 
quadratic  regulator  performance  measure  is  modified  to  allow 
for  its  application  to  the  tracking  problem  solution. 
Additionally,  the  control  systems  are  designed  using  reduced 
order  models.  Computer  simulations  show  that  the  reduced 
order  controller  designs  provide  results  comparable  to  the 
full  order  controller  designs.  Successful  hardware  tests  with 
roll  rate  control  system  validated  the  reduced  order  model 
design  philosophy  used  in  this  research. 
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I .  INTRODUCTION 


A.  THE  ARCHYTAS  CONCEPT 

The  current  inventory  of  Unmanned  Air  Vehicles  (UAVs)  is 
not  able  to  meet  the  expanding  need  for  real-time  intelligence 
at  the  Marine  expeditionary  unit  or  Army  battalion  level.  The 
Naval  Postgraduate  School  UAV  Flight  Research  Lab  is  directing 
efforts  at  developing  a  ducted-fan  vertical  takeoff  and 
landing  (VTOL)  vehicle  to  meet  these  increasing  needs.  The 
NPS  air  vehicle,  named  Archytas,  is  serving  as  a  technology 
demonstrator  to  evaluate  the  concept  of  a  winged  ducted-fan 
VTOL  aircraft.  The  research  is  being  directed  at  applying  the 
technology  and  equipment  developed  in  the  U.S.  Marine  Corps' 
Airborne  Remotely  Operated  Device  (AROD)  program  and  the  U.S. 
Army's  AQUILA  program. 

Archytas,  pictured  in  Figure  1.1,  is  designed  to  take  off 
and  land  vertically.  After  climbing  to  altitude,  Archytas 
will  transition  to  horizontal  flight  by  pitching  about  its 
center  of  gravity  to  a  wings  level  attitude.  The  positioning 
of  the  duct  and  wings  ( including  the  canard)  allow  for  the 
vertical  takeoff  and  landing  capability.  The  ability  to 
transition  to  horizontal  flight  will  extend  the  vehicle's 
range  and  provide  the  capability  for  a  high  speed  dash. 
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TOP  VIEW 


Figuro  1.1  Sketch  of  Arehytaa 


B.  THE  CONTROL  PROBLEM 

Archytas  is  powered  by  a  vertically  mounted  28-horsepower 
engine  turning  a  three-bladed  propeller.  The  use  of  a  single 
propeller  in  a  duct  (ducted  fan)  simplifies  the  design,  but 
intensifies  the  stability  and  control  problems.  The  dynamic 
behavior  about  a  given  axis  is  coupled  with  other  vehicle 
dynamics,  in  particular,  there  are  three  types  of  coupling: 

1.  The  single  propeller  design  introduces  a  gyroscopic 
coupling  between  the  pitch  and  yaw  axes. 

2.  Reactive  torques  are  applied  to  the  roll  axis  as  the 
engine  speed  is  varied. 

3.  A  loss  of  lift  due  to  thrust  occurs  when  the  vehicle  is 
pitched  during  the  translation  to  horizontal  flight. 

It  was  the  goal  of  this  thesis  to  design  a  control  system 
which  would  allow  stable  flight  of  Archytas  during  takeoff  and 
landing,  and  during  the  vertical  mode  of  operation.  Because 
of  the  coupled  nature  of  the  Archytas  control  problem,  linear 
quadratic  regulator  control  theory  was  used. 

c.  THESIS  ORGANIZATION 

In  Chapter  II,  the  general  nonlinear  equations  of  motion 
are  developed.  These  equations  of  motion  are  then  applied  to 
Archytas  with  special  attention  given  to  the  effects  of  the 
spinning  propeller.  Additionally,  the  equations  of  motion  are 
linearized  about  the  hover  operating  condition  using  the 
small-disturbance  theory. 
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Chapter  III  discusses  the  state  space  representation  of  a 
general  system  and  develops  a  procedure  for  selecting  state 
variables  for  the  optimal  control  tracking  problem.  Key  to 
the  design  of  the  Archytas  control  system  is  the  proper 
selection  of  the  state  variables. 

In  Chapter  IV,  three  control  law  designs  are  formulated 
based  on  the  linearized  hover  equations.  One  design  is  for 
the  single  input,  single  output  (SISO)  roll  rate  controller. 
The  second  design  is  for  the  SISO  altitude  rate  controller. 
These  two  SISO  systems  are  similar  in  their  development. 
Next,  the  multiple  input,  multiple  output  (MIMO)  pitch  and  yaw 
angle  controller  is  designed.  Central  to  each  control  law 
design  is  the  use  of  a  reduced  order  model  to  simplify  the 
design  process  and  physical  implementation.  Finally,  these 
control  laws  are  applied  to  the  nonlinear  model  for 
validation. 

Chapter  V  discusses  the  field  test  results  of  the  roll 
rate  controller.  The  roll  rate  controller  was  evaluated  with 
Archytas  mounted  on  a  test  stand  to  allow  a  roll  about  the 
longitudinal  axis.  In  addition,  conclusions  based  on  the 
computer  simulations,  the  field  tests,  and  recommendations  for 
future  research  are  presented. 
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II.  MODELING  THE  ARCHYTAS  PROTOTYPE 

The  purpose  of  this  chapter  is  to  develop  a  suitable 
dynamic  model  of  the  Archytas  prototype.  Because  Archytas 
is  a  ducted  fan  device,  special  attention  must  be  given  to 
the  significant  gyroscopic  contribution  of  its  propeller. 

A.  DERIVATION  OP  RIGID  BODY  EQUATIONS  OF  MOTION 

The  rigid  body  equations  of  motion  in  this  section  are 
developed  for  the  Archytas  prototype  in  the  following  way. 
First,  Archytas  is  regarded  as  a  single  rigid  body,  and  the 
equations  of  motion  are  derived  with  respect  to  a  set  of 
body  fixed  axes.  These  equations  are  the  general  equations 
governing  aerodynamic  flight  for  all  aircraft.  Next,  the 
changes  introduced  by  the  spinning  rotor  are  evaluated  and 
included  in  the  equations.  These  are  the  gyroscopic  effects 
due  to  the  propeller.  Finally,  the  development  of  a 
complete  model  is  undertaken  for  the  specific  case  of 
Archytas  using  the  actual  measurements  and  experimental  data 
from  the  AROD  prototype  as  first  approximations. 

l*  Force  and  Moment  Equations 

The  general  equations  of  motion  are  developed  for  a 
typical  aircraft  in  References  l  and  2.  A  combination  of 
the  two  approaches  is  taken  here  to  arrive  at  the  set  of 
equations  describing  Archytas.  The  equations  of  motion  are 
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obtained  from  Newton's  second  law,  which  states:  The 
summation  of  all  external  forces  acting  on  a  body  is  equal 
to  the  time  rate  of  change  of  the  momentum  of  the  body;  and 
the  summation  of  the  external  moments  acting  on  the  body  is 
equal  to  the  time  rate  of  change  of  the  moment  of  the 
momentum  (angular  momentum) .  The  time  rates  of  change  of 
linear  and  angular  momentum  are  referred  to  an  absolute  or 
inertial  reference  frame.  This  absolute  or  inertial 
reference  frame  is  an  axis  system  fixed  to  the  Earth. 

Figure  2.1  depicts  both  the  body  fixed  axes  and  the  inertial 
reference  frame. [Ref.  1:  p.84] 

Newton's  second  law  can  be  expressed  in  the 
following  vector  equations: 

;  (2.1) 

<2-2> 

where  £  is  the  externally  applied  force,  M  the  externally 
applied  moment  about  the  center  of  mass,  v  the  velocity 
vector,  and  H  the  angular  momentum  vector  about  the  center 
of  mass. 

The  vector  equations,  in  scalar  form,  consist  of 
three  force  equations  and  three  moment  equations.  The  force 
equations  can  be  expressed  as 
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(2.3) 


>  Fy'-§i{mv)  '  1 

where  Fx,  Fyl  Fz  and  u,  v,  w  are  the  components  of  the  force 
and  velocity  along  the  x,  y  and  z  axes,  respectively.  The 
force  components  are  composed  of  contributions  due  to  the 
aerodynamic,  propulsive,  and  gravitational  forces  acting  on 
the  aircraft,  the  Archytas  prototype  for  the  purpose  of  this 
thesis.  The  moment  equations  can  be  expressed  in  a  similar 
manner : 

;  ”-&*•’  (2-4’ 

where  L,  M,  N  and  Hy,  Hz  are  the  components  of  the  moment 
and  moment  of  momentum  along  the  x,  y  and  z  axes 
respectively. 

Now,  considering  Figure  2.2,  let  5m  be  an  element  of 
mass  of  the  Archytas  prototype,  v  be  the  velocity  of  the 
elemental  mass  relative  to  the  inertial  axes,  and  let  5£  be 
the  resultant  force  that  acts  upon  it.  Newton's  second  law 
then  gives  the  equation  of  motion  of  5m: 

.  (2-S) 

at 

The  total  force  acting  on  the  vehicle  is  a  summation  of  all 
the  forces  that  act  upon  all  the  elements.  The  internal 
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forces,  those  exerted  by  one  element  upon  another,  all  occur 
in  equal  and  opposite  pairs,  by  Newton's  third  law,  and 
hence  contribute  nothing  to  the  summation.  Thus  L6 £=£  is 
the  resultant  external  force  acting  upon  the  vehicle.  The 
velocity  of  the  differential  mass  6m  is  given  by 

f  <*•«> 

c  dt 

were  v,.  is  the  velocity  of  the  center  of  mass  of  the 
aircraft  and  d£/dt  is  the  velocity  of  the  element  relative 
to  the  center  of  mass.  Substituting  this  expression  for  the 
velocity  into  Newton's  second  law,  equation  (2.1),  yields 

E6£=£=£5jn-jL(*c+i|)  .  (2.7) 

Assuming  that  the  mass  of  the  aircraft  is  constant,  equation 
(2.7)  can  be  rewritten  as 


z=4 

dt 


E  <*c 


♦I’6” 


(2.8) 


or 


_  dy 
E=zn  c 


dt  dt 


(2.9) 


Because  £  is  measured  from  the  center  of  mass,  the  summation 
li.6m  is  equal  to  zero  and  the  force  equation  (2.9)  becomes 
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(2.10) 


z~  ~dt  ' 

which  relates  the  external  force  on  the  aircraft  to  the 
motion  of  the  vehicle's  center  of  mass. 

The  relationship  between  the  external  moment  and  the 
rotation  of  the  aircraft  is  obtained  from  a  consideration  of 
the  moment  of  momentum.  For  the  differential  element  of 
mass,  6m,  the  moment  of  momentum  is  by  definition  6ff=j:xv6;n. 


The  moment  equation  can  be  written  as 


bM=-£z&H=-£-  (£XJf)  S/n  . 
dt  dt 


(2.11) 


The  velocity  of  the  mass  element  can  be  expressed  in  terms 
of  the  velocity  of  the  center  of  mass  and  the  velocity  of 
the  mass  element  relative  to  the  center  of  mass: 


V=V  +j£X £  ; 

x  dt~c 


(2.12) 


where  a  is  the  angular  velocity  of  the  vehicle  and  r  is  the 
position  of  the  mass  element  measured  from  the  center  of 
mass.  The  total  moment  of  momentum  can  be  written  as 

5jn+E  (Siixz)  ]  5m  .  (2.13) 


The  velocity  ^  is  a  constant  with  respect  to  the  summation 
and  can  be  taken  outside  of  the  summation  sign 
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(2.14) 


The  first  term  in  equation  (2.14)  is  zero  because  the  term 
Ut6m=0  as  explained  previously.  The  position  vectors  and 

angular  velocity  can  be  expressed  as 

£=*jL+yj+z&  ;  (2.i5) 

id.=p±+qj+r&  ;  (2.16) 

where  p,  q  and  r  are  the  scalar  components  of  and  A,  j,  )£ 
are  unit  vectors  in  the  directions  of  x,  y  and  z. 
Substituting  a  and  r  into  equation  (2.14)  and  expanding,  H 
can  be  written  as 

(Pl+ql+rjk)  £  (x2+y2+z2)5jn-]T  (xA+yi+z&)  (px+qy+rz)  t>m  . 


The  scalar  components  of  H  are 


(2.17) 


Hx=  p£(y2+z2)  Sm-qYjXy  5 m-r^xz  tun  ; 

Hy=-pYjXy  5/n+g£(x2  +  z2)  tun  -zYyz  & m  !  (2.18) 

Hz=-pY,xz  tun-q^yz  5r»+r£(x2+y2)  5/n  . 

The  summations  in  the  above  equations  are  the  mass  moment 
and  products  of  inertia  of  the  aircraft  and  are  defined  as 


12 


(2.19) 


-V-///  (y2  +  Z2)5m  ,  ~xy~ff f*y  5 m  ; 

V  V 

Iy=fjf  (x2  +  z2y>*m  '  Zxz=SifXZ  5jn  ; 

V  V 

-,■///  (x2+y2)  jfyz  • 

V  V 

The  terms  Ix,  iy  and  Iz  are  the  mass  moments  of  inertia  of 
the  body  about  the  x,  y  and  z  axes,  respectively.  The  terms 
with  the  mixed  indices  are  called  the  products  of  inertia. 
Both  the  moments  and  products  of  inertia  depend  on  the  shape 
of  the  body  and  the  manner  in  which  its  mass  is  distributed. 
The  larger  the  moments  of  inertia  the  greater  the  resistance 
the  body  will  have  to  rotation.  Applying  the  notation  of 
equation  (2.19)  to  equation  (2.18),  yields  the  scalar 
equations  for  the  moment  of  momentum: 

^x~  P^x~^Zxy~z^xz  • 

Hy  =-pIxy+ qly - r Iyx  ;  (2.20) 

pi xx  ~  Qlyz + Z I x  * 

If  the  reference  frame  is  not  fixed  to  the  aircraft,  then, 
as  the  aircraft  rotates,  the  moments  and  products  of  inertia 
will  vary  with  time.  To  avoid  this  difficulty  an  axis 
system  will  be  fixed  to  the  aircraft  (body  axis  system) . 

Now  the  derivatives  of  the  vectors  v  and  H  referred  to  the 
rotating  body  frame  of  reference  must  be  determined. 
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It  can  be  shown  that  the  derivatives  of  an  arbitrary 
vector  A  referred  to  a  rotating  body  frame  having  an  angular 
velocity  can  be  represented  by  the  following  vector 
identity 


d& 

dt  z 


d& 

dt  B 


(2.21) 


where  the  subscripts  I  and  B  refer  to  the  inertial  and  body 
fixed  frames  of  reference,  respectively.  Applying  this 
identity  to  equations  (2.1)  and  (2.2)  yields 


d% 

dt  n  c 


(2.22) 


.  (2.23) 

These  are  the  general  equations  governing  aerodynamic  flight 
and  have  the  scalar  components: 

Fx=m{u+qw-iv)  ,  Fy=m(-v+iu-pw)  ,  Fz=m(W+pv-qu)  ;  (2.24) 

L~Hx+qHz-rHy,  M=Hy + rHx ~pHt ,  N=Hz+pHy~qHx  .  (2.25) 

The  components  of  the  force  and  moment  acting  on  the 
aircraft  are  composed  of  aerodynamic,  gravitational  and 
propulsive  contributions. 

At  this  point,  it  is  recognized  that  most  aircraft 
have  a  plane  of  symmetry.  If  the  xz  plane  is  selected  to 
coincide  with  this  plane  of  symmetry,  then  from  equation 
(2.19),  I ^=1^=0  must  be  satisfied.  However,  for  the  case  of 
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the  Archytas,  the  three  blades  of  the  propeller  provide  two 
planes  of  symmetry,  xz  and  xy.  Thus,  the  products  of 
inertia,  1^,  Iyz  and  I^,  equal  zero.  The  moment  equations 
for  Archytas  can  now  be  written  as 

L=IJ>*qr  { Ig-Iy)  ; 

Af=  Iy<Z+rp ( Jx- J2 )  ;  (2.26) 

N=Izt+pq(Iy-Ix)  . 

2.  Effect  of  the  Spinning  Rotor 

Archytas,  like  AROD,  is  a  gyroscope.  The  single 
propeller  rotates  about  the  longitudinal  vehicle  axis  to 
produce  a  downwash  or  jet  of  air  through  the  duct  which 
makes  up  the  Archytas  body.  This  spinning  rotor  exerts  a 
gyroscopic  moment  on  the  body  of  the  vehicle.  Reference  3 
states  that  in  developing  the  equations  of  motion  for 
aircraft  with  propellers  which  exert  gyroscopic  moments  on 
the  body  "more  often  than  not,  such  gyroscopic  moments  turn 
out  to  be  negligible."  However,  as  demonstrated  by  Bassett 
[Ref.  4:  p.  19],  in  the  AROD  case  the  angular  momentum  of 
the  propeller,  Ip(AROdh  equals  11.3  ft2-lbm/sec.  Compared 
with  AROD' s  nominal  total  mass  of  2.64  lbm  (85  lbw)  ,  it  is 
clear  that  the  angular  momentum  imparted  by  the  propeller  is 
significant  and  that  gyroscopic  effects  will  play  a  large 
part  in  modeling  the  dynamic  behavior  of  AROD.  For 
Archytas,  Ip  is  identically  equal  to  Ip(AROD)f  11  • 3  ft2-lbm/sec. 
Compared  with  Archytas'  nominal  total  mass  of  3.11  lbm  (100 


lbw) ,  it  is  clear  that  similar  to  AROD  the  gyroscopic 
effects  will  be  significant  and  must  be  included  in  modeling 
the  dynamic  behavior  of  Archytas.  This  gyroscopic  moment 
can  be  accounted  for  as  follows. 

Angular  momentum,  Hp,  due  to  the  propeller  is 
defined  as 

Hp-Ip<*p=lHpx+jHpy+l&pz  ;  (2.27) 

where  lp  is  the  propeller  moment  of  inertia,  and  a>p  is  the 
propeller  angular  velocity,  since  the  propeller  lies  in  the 
yz  plane  and  spins  symmetrically  about  the  x  axis,  is 

directed  only  along  x  and  H^H^O.  Equation  (2.27)  becomes 

Hp=IpUp=iHpx  .  (2.28) 

Etkin  [Ref.  2:  p.93]  states  that  the  resultant  angular 
momentum  of  an  aircraft  with  spinning  propellers  is  obtained 
simply  by  adding  H,,  to  the  IJ  previously  defined  by  equation 
(2.20).  Adding  equations  (2.28)  and  (2.20)  and  keeping  in 
mind  that  the  products  of  inertia  equal  zero,  yields 

Hx=plx+  IpWp  ; 

Hy=qly  ;  (2.29) 

Hz=rlz  . 

Applying  equation  (2.29)  to  equation  (2.25),  the  moment 
equations  can  be  written  for  the  specific  case  of  Archytas 
as 
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L=IJ>+qr(I3-Iy)  +Ipu>p  ; 

Af=  Iyg*  rp{Ix~Iz)  +  r  JpQp  ;  (2.30) 

N=Izr+pq(Iy-Ix)  -gJp«p  ; 

and  the  force  equations  are  those  of  equation  (2.24). 

3.  Orientation  and  Position 

Because  the  frame  of  reference  developed  for  the 
equations  of  motion  is  fixed  to  the  aircraft,  and  moves  with 
it,  the  position  of  the  aircraft  cannot  be  described 
relative  to  it.  The  orientation  and  position  of  the 
aircraft  can  be  defined  in  terms  of  a  fixed  frame  of 
reference  as  shown  in  Figure  2.3.  [Ref.  1:  p.  89] 

The  orientation  of  the  aircraft  can  be  described  by 
a  series  of  three  consecutive  rotations,  whose  order  is 
important.  The  angular  rotations  are  called  the  Euler 
angles.  The  orientation  of  the  body  frame  with  respect  to 
the  fixed  frame  can  be  determined  in  the  following  manner. 
The  aircraft  is  imagined  first  to  be  oriented  so  that  its 
axes  are  parallel  to  the  fixed  frame  and  the  following 
rotations  are  then  applied.  [Ref.  2:  pp.  89-91] 

1.  A  rotation  ¥  about  oz,,  carrying  the  axes  to  Cx2y2z2 
(bringing  Cx  to  its  final  azimuth) . 

2.  A  rotation  0  about  oy2,  carrying  the  axes  to  Cx3y3z3 
(bringing  Cx  to  its  final  elevation) . 

3.  A  rotation  $  about  ox3,  carrying  the  axes  to  their 
final  position  Cxyz  (giving  the  final  angle  of  bank  to 
the  wings) . 
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Now  that  the  Euler  angles  are  defined,  the  flight  velocity 
components  relative  to  the  earth-fixed  reference  frame  can 
be  determined.  To  do  this,  let  the  velocity  components 
along  the  x,,  y,,  z,  frame  be  dx/dt,  dy/dt,  dz/dt  and, 
similarly,  let  the  subscripts  2  and  3  denote  the  components 
along  x2,  y2,  z2  and  x3,  y3,  z3,  respectively.  From  Figure 
2.3,  it  can  be  shown  that 


where 


and 


and 


dx 

— —  =u, 
dt  1 


dt  1 


dz 

dt 


=w. 


(2.31) 


14=05  cos  tlr-v2sini|r  ; 

vx=u2sinijf+v2cosT|f  ;  (2.32) 


u2=u3cos0+w3sin0  ; 

V2  =  V3  ;  (2.33) 

w2  =  -OjSin0  +  w3cos0  ; 

u3  =  u  ; 

v3=vx:os®-sin®  ;  (2.34) 

w3  =  vsin<&+wcos<t  ; 


from  this,  the  absolute  velocity  in  terms  of  the  Euler 
angles  and  velocity  components  in  the  body  frame  can  be 
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determined.  Note  the  shorthand  notation  S^ssin^,  c^scosS^, 
Sessin0,  etc,  used  in  the  following  equations: 


dx 

dt 

dy 

dt 

dz 

cftj 

. 

CeC^ 

-St 


5$  Cy 


+SqS^ 

u 

-SqC^ 

V 

C+Cq 

w 

(2.35) 


Integration  of  these  equations  yields  the  aircraft's 
position  relative  to  the  fixed  frame  of  reference.  The 
relationship  between  the  angular  velocities  in  the  body 
frame  (p,  q  and  r)  and  the  Euler  rates  (  0 ,  i|»  and  )  can 

also  be  determined  from  Figure  2.3. 


p 

■1  0  -se 

4 

Q 

= 

0  c* 

6 

r 

0  -Sq  cecy 

* 

(2.36) 


Equations  (2.36)  can  be  solved  for  the  Euler  rates  in  terms 
of  the  body  angular  velocities  and  is  given  by  equation 
(2,37) 


By  integrating 
and  $)  can  be 


4 

6 

= 

1  s*tan0  C*tan0 
0  C*  ~S, 

0  S*sec0  C*sec0 


P 

Q 

r 

(2.37) 


the  above  equations,  the  Euler  angles  (0, 
determined. 


4.  Gravitational  and  Thrust  Forcas 

The  gravitational  force  acting  on  the  aircraft  acts 
through  the  center  of  gravity  of  the  aircraft.  Because  the 
body  axis  system  is  fixed  to  the  center  of  gravity,  the 
gravitational  force  will  not  produce  any  moments.  However, 
the  gravitational  force  will  contribute  to  the  external 
force  acting  on  the  aircraft  and  will  have  components  along 
the  respective  body  axes.  From  Figure  2.4  the  gravitational 
force  component  in  the  direction  of  each  axis  is  found  to  be 

Xg=-mg  cosflcosY  ; 

Yg=  mg  cosOsin’P  ;  (2.38) 

Zg=  sinfl  . 

With  the  aerodynamic  forces  (including  the  propulsive 
forces)  denoted  by  (X,  Y,  Z) ,  the  resultant  external  forces 
are 

‘  Fx=X-mg  cosOcos’P ; 

Fy=Y+lhg  cosflsinT;  (2.39) 

Fz=Z-mg  sin0  . 

5.  summary  of  Equations  of  Notion 

In  the  previous  sections,  the  equations  that 
completely  describe  the  dynamic  behavior  of  Archytas  have 
been  developed.  Equations  (2.24)  and  (2.30)  define  the 
externally  applied  forces  and  moments  which  are  represented 
by  Fx,  Fy ,  Fx  and  L,  M,  N.  Through  the  Euler  angles,  the 
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behavior  of  Archytas  can  be  observed  relative  to  the  Earth. 
Specifically,  the  translational  velocities  for  the  fixed 
frame  of  reference,  dx/dt,  dy/dt  and  dz/dt,  can  be 
determined  from  the  body-fixed  velocities,  u,v,w  and  the 
Euler  angles,  0,  ¥  and  t,  using  the  transformation  of 
equation  (2.35).  Additionally,  equation  (2.37)  describes 
the  relationship  between  the  Euler  angles  and  the  body 
angular  velocities,  p,  q  and  r.  Table  2.1  gives  a  summary 
of  the  rigid  body  equations  of  motion. 

B.  ARCHYTAS  NONLINEAR  SYSTEM  EQUATIONS 
1.  Applied  Forces  and  Moments 

The  Archytas  model  is  developed  using  the 
measurements  and  data  from  the  AROD  prototype  as  first 
approximations.  These  measurements  and  data  were  obtained 
by  the  AROD  project  engineers  based  on  wind  tunnel  tests. 
This  data  consists  of  tabular  results  that  describe  the 
aerodynamic  lift  and  drag  coefficients  and  physical 
measurements  of  constants  such  as  weight,  moments  of  inertia 
and  servo  gains.  This  tabulated  data  forms  the  basis  from 
which  the  applied  forces  and  moments  may  be  determined.  The 
forces  and  moments,  which  are  computed  from  the  data  are  of 
two  types:  aerodynamic  and  thrust.  This  data  is  listed  in 
Appendix  A. 
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TABLE  2.1  SUMMARY  OP  GENERAL  EQUATIONS  OF  MOTION 


X-mg  C9Cy  =  m(u+qw-iv) 
Y+mg  C9Sy  =  mi  v+ru-pw) 
Z-mg  S6  -  m{ti+pv~qu) 

L~Ij>*qz(Ix-Iy)*Ipi»I> 
M~Iyq*Tp{Ix-Ix)  +rlpap 
N-Izi+pq(ly-lx)  -qlpu>p 

p«#-TSe 
<7=eC*+'PC95» 
r=*C8C,-dS, 

0=gC*-r5* 
<fr-p+gs*r8+rc,r8 
't=(qS9*rC9)  seed 

Velocity  of  aircraft 
angles  and  body 

dx 1 
dt! 
dy 
dt 
dz 
dt 


Force  equations 


Moment  equations 


Angular  velocities 


Euler  rates 


in  the  fixed  frame  in  terms  of  Euler 
velocity  components 


C8S8C8+S#S8 

-58  5jC8  C*C8 


a.  Total  Angle  of  Attack  and  Body  Roll  Angle 

The  aerodynamic  data  describes  the  forces  and 
moments  relative  to  the  vehicle's  total  velocity  vector, 
vTOT.  These  forces  and  moments  (in  the  body-fixed  coordinate 
system)  depend  on  the  total  angle  of  attack  (a')  and  the 
body  roll  angle  (A) .  The  total  angle  of  attack  and  body 
roll  angle  can  be  defined  in  terms  of  the  velocity 


components  as  shown  in  Figure  2.5.  The  equations  for  a '  and 
A  are  given  as: 


and 


a -=sin -.««! 

u 


(2.40) 


A=tan_1—  .  (2.41) 

w 

b.  Aerodynamic  Forces  and  Moments 

The  forces  computed  from  the  tabular  data  are 
lift  (F,)  and  drag  (Fd) .  These  forces  are  depicted  in  Figure 
2.6(a).  The  transformation  from  lift  and  drag  to  Fu,  F,y  and 
F„  is  given  as 

*Vx  _  sina'  -cosa'  Fi  _  (2.42) 

Fayz  [-cosa  '  -sina  'J  Fa\ 

and 


Fay  fsinAl  (2.43) 

Faz\  ‘  [cosAj  No¬ 
where  F„,  Fty  and  Fu  are  the  forces  in  the  x,  y  and  z 
directions  due  to  the  aerodynamic  data. 

Similarly,  the  aerodynamic  moments  applied  to  the 
body  axes  as  a  result  of  the  vehicle7 s  movement  through  the 
air  can  be  derived.  These  moments  (shown  in  Figure  2.6(b)) 
are  referred  to  as  the  aerodynamic  angular  moments  of  roll 
(LJ  ,  pitch  (MJ  and  yaw  (NJ  and  are  given  as 
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y„-zb  plane 


z 
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(2.44) 


’  La 

= 

sina  ' 

-cosa  '] 

My 

Mpyyz 

-cosa  ' 

-sina  'j 

Mr. 

and 


Ma 

sinA 

cosA 

where  My,  and  M,  are  the  yaw  and  roll  moments  relative  to 
vtot*  is  the  moment  in  the  y-z  plane.  The  relationships 

that  define  F,,  Fd,  and  My,  M,  were  developed  by  the  AROD 
project  engineers  and  are  listed  in  Appendix  A.  [Ref.  15] 

c.  Control  Forces  and  Moments  Due  to  Command  Inputs 

The  forces  and  moments  previously  discussed  are  a 
result  of  the  vehicles  motion  through  the  air.  A  second 
category  includes  those  forces  and  moments  which  are  a 
result  of  the  commanded  inputs.  These  inputs  control  the 
rotor  speed  and  the  displacement  of  the  control  surfaces. 
[Ref.  4:  p.  33] 

(1)  Forces  Due  to  Induced  Thrust  of  the  Ducted 
Fan.  The  commanded  inputs  include  the  ability  to  change  the 
rotor  speed;  thus,  changing  the  thrust  provided  by  the 
ducted  fan.  Due  to  the  orientation  of  the  body-fixed  axes, 
the  force  due  to  the  thrust  (F^)  is  directed  completely 
along  the  x-axis.  The  relationship  that  defines  F^  was 
developed  by  the  AROD  project  engineers  and  is  given  as 
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Frhr  -  THOVER  +  XRPM  *  Sr?71  (2.46) 

where  THOVER  is  a  constant  equal  to  the  nominal  thrust  at 
hover.  THOVER  is  set  equal  to  Archytas  prototype  weight  of 
76.5  ft/lbs.  XRPM  is  the  slope  of  the  thrust  versus  engine 
rpm  curve.  For  the  hover  rpm  of  712.0943  rad/sec  (6800 
rpm)  ,  XRPM  is  equal  to  0.2387  lbf/rad/sec.  bzpm  is  the 

change  in  engine  rpm  from  the  nominal  hover  rpm  in  rad/ sec. 

(2)  Moment  Due  to  Ducted  Fan  Effects .  A 
gyroscope  imparts  no  torque  on  its  axis  if  it  spins  with  a 
constant  angular  rotation.  In  hover  (constant  rotor  speed) , 
Archytas  behaves  similarly  to  a  gyroscope.  If  the  rotor 
accelerates  (positively  or  negatively)  a  torque  is  applied 
to  the  axis.  This  torque  is  accounted  for  in  Equation 
(2.30).  However,  Archytas  is  a  ducted  fan  and  the  drag 
between  the  rotor  tip  and  the  inside  body  wall  creates  a 
moment  about  the  x-axis  (roll) .  The  project  engineers  for 
AROD  determined  an  approximation  for  this  moment  based  on 
experimentation.  Because  the  Archytas  duct  is  identical  to 
the  AROD  duct,  the  moment  determined  for  AROD  applies  to 
Archytas  and  is  given  as 

Thr  =  Kductbrpm  (2.47) 

where  6rpm  is  defined  above.  is  a  constant  which  is 

dependent  on  the  duct  geometry  and  is  equal  to  0.072Q.  [Ref. 
4:  pp.  33-34] 
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(3)  Moments  Due  to  Control  Surface  Displacements . 
The  commanded  inputs  also  include  the  displacement  of  the 
four  control  surfaces  (vanes)  within  the  downwash  from  the 
duct.  The  displacement  of  these  vanes  within  the  downwash 
of  the  duct  imparts  moments  about  the  body  axes.  Figure  2.7 
shows  how  the  vanes  are  symmetrically  arranged  below  the 
duct.  The  vanes  are  displaced  by  a  servo  mechanism 
connected  directly  to  the  top  of  each  vane.  Vanes  (1)  and 
(3)  are  operated  together  as  "elevators'*  and  impart  a  moment 
about  the  y-axis  (pitch).  Vanes  (2)  and  (4)  together  are 
the  "rudder”  and  contribute  a  moment  about  the  z-axis  (yaw) . 
Vanes  (1)  and  (3)  displaced  in  opposite  directions  and  (2) 
and  (4)  displaced  oppositely  work  as  "ailerons"  to  impart  a 
moment  about  the  x-axis  (roll) .  The  actual  torque  applied 
by  each  combination  of  vanes  was  determined  experimentally 
and  described  by  "constants  of  effectiveness”  which  were 
calculated  by  the  AROD  project  engineers.  These  constants 
of  effectiveness  can  be  applied  directly  to  Archytas.  [Ref. 

4:  p.34] 

The  constants  of  effectiveness  are  given  the  symbols 
L»ff»  Mjef,  and  Nreff,  for  their  contribution  of  moments  about  the 
roll,  pitch  and  yaw  axes  due  to  a  displacement  by  the 
ailerons,  elevator  and  rudder.  The  relationships  resulting 
in  moments  about  the  three  body  axes  are 
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Figure  2.7  Archytas  control  Van as 
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(2.48) 


^ faccoz^  e  ' 

^C~^reff^  factor^  z  ' 

where  ba,  6g  and  6r  are  the  displacements  of  the  aileron, 

elevator  and  rudder  respectively.  Pftctor  and  YfKUx  are  scaling 
factors.  The  relationships  that  define  Pflctt)r  and  Yfaclor  were 
developed  by  the  AROD  project  engineers  and  are  listed  in 
Appendix  A.  L^,  M«ff  and  Nreff  equal  -150,379.57,  -112,716.87 
and  -128,774.80,  respectively,  with  units  of  lbm-in2- 
rad/sec2.  Because  of  the  symmetry  of  Archytas,  cross 
coupling  of  the  control  surfaces  is  negligible  and  is 
ignored.  [Ref.  4:  pp.  35-36] 

2.  servo  Equations  for  Control  Surfaces  and  Throttle 
The  model  airplane  servos  used  to  drive  the  control 
vanes  and  throttle  linkage  are  identical  and  can  be  modeled 
as  second-order  dynamical  systems.  The  response  of  these 
servos  to  a  step  response  was  measured.  These  measurements 
were  used  to  compute  the  natural  frequency  and  damping 
ratio.  The  results  are  summarized  below  in  the  form  of 
constants  H,  and  H2.  A  detailed  explanation  as  to  how  these 
results  were  obtained  is  contained  in  Appendix  B. 
a.  Control  Surface  Servos 

Three  equations  will  describe  the  operation  of 
the  elevators,  rudders  and  ailerons.  For  each  of  these 
equations,  at  least  two  servos  are  operating  at  the  same 
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time  on  different  control  vanes.  Each  servo  receives  a 
command  input  (pulse-width  modulated  signal)  which  results 
in  an  angular  displacement  of  the  servo.  The  corresponding 
differential  equations  for  the  servos  can  be  written  as 

Sa  =  -tf6a-//25a+tf2ua  ; 

hB=-H.tg-H2h9+H2U9  ;  (2.4‘ 

<!>  Z=-H1& T~H2t> z+H2Ur  ; 

where  H,  and  H,  are  the  servo  gain  constants  equal  to  71.1 
and  2745.8.  u,,  ue  and  ur  are  the  servo  inputs.  6a,  6*  and 

5r  are  the  servo  position  angles. 

b.  Throttle  Servos  and  Engine 

The  servo  motor  used  to  open  and  close  the 
throttle  is  identical  to  the  servo  motors  used  for  the 
control  surfaces.  The  2-cycle,  2-cylinder  gasoline  engine 
with  dual  carburetors  can  be  modeled  as  a  first  order  lag 
system.  The  complete  third  order  system  is 


5t=-H14t.-fl2d{:+tf2ut.  ; 

^rpm=_0)B^rpm  +  <,)£^E^r  ' 


(2.50) 


where  is  the  lag  time  constant  equal  to  2.0  rad/ sec  and 
Ke  is  a  scaling  factor  equal  to  837.758  rad/sec/rad. 


3 .  Summary 

The  result  of  the  previous  section  was  nine 
equations  completely  describing  the  dynamic  behavior  of 
Archytas.  These  equations,  combined  with  the  equations  that 
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define  the  servos  and  engine,  form  the  complete  nonlinear 
system.  This  nonlinear  system  will  be  the  base  model 
throughout  this  thesis.  Table  2.2  lists  the  equations 
rearranged  so  that  the  dynamic  variables  of  interest  may  be 
determined.  Additionally,  the  applied  forces  and  moments 
defined  in  this  section  have  been  substituted  into  the 
equations.  These  applied  forces  and  moments  complete  the 
development  of  the  model  for  the  specific  case  of  Archytas. 

c.  LINEARIZING  THE  ARCHYTAS  MODEL 

There  exist  many  analytical  and  graphical  techniques  for 
controller  design  and  analysis  of  linear  systems. 

Conversely,  there  are  no  good  methods  available  for  solving 
a  wide  class  of  nonlinear  systems.  Thus,  in  the  design  of 
control  systems  it  is  practical  to  first  design  the 
controller  based  on  the  linear  system  model  generated  by 
neglecting  the  nonlinearities  of  the  system.  These 
nonlinearities  are  neglected  by  linearizing  the  model  about 
a  steady-state  reference  condition.  The  designed  controller 
is  then  applied  to  the  nonlinear  system  model  for  validation 
and  subsequent  redesign  if  necessary.  In  this  section,  a 
linear  model  is  generated  from  the  equations  of  Table  2.2 
based  on  steady-state  assumptions  and  physical 
approximations.  [Ref.  5:  p.  nj 

The  nonlinearities  of  the  Archytas  system  equations  fall 
into  two  categories:  (l)  nonlinear  combination  of  states 
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_ Table  2.2  ARCHYTA8  NONLINEAR  SYSTEM  EQUATIONS _ 

J>u4-  Ua(MylMr)  +LTht(6  )  +Lt(Sa)  +(ly-lz)rq-ITxto  ] 

<7=4*  [Ma{My,M:)  +Afc(5&)  +{IS-Ix)pr-Izxupr] 

1y 

i~~j~  WMylMt)  +JVt(5r)  +(ix-iy)pq+itxupq) 

±Z 

u=  (rv-gv)  -gcosecosV  +  A  iF^tF^F*)  +FThrx{ 6rpm) ) 

v=  (pv-ru)  +gcos0sin’F  +  iF.v(F,,  Fd) 

m  y 

&={qu-pv)  -grsinQ + -i Faz ( F1 ,  Fd) 

8=qcos<&-rsin4> 

<t  =p + grs  i  n®  t  an0 + r  c  o  s  ®  t  and 
¥  =  (gsin®  +rcos$)  sec0 

8  r  ^.8  r-^2®r+^2Ur 
8 t— «i8  t-H2*t+H2ut 
8  rpm_~ws8rpm+wff^jr8t 

(e.g.  Equation  (2.30))  and  (2)  discontinuous  functions  (e.g. 
table  lookup  of  aerodynamic  force  and  moment  coefficients) . 
These  nonlinearities  will  be  neglected  and  the  nonlinear 
equations  will  be  replaced  with  linear  approximations. 

1.  Steady-State  Assumptions 

The  motion  of  Archytas  in  the  hover  mode  consists  of 
small  perturbations  from  a  steady-state  condition.  The 
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steady-state  hover  condition  is  defined  as  a '  equal  to  90°. 
All  translational  and  angular  movement  is  very  small.  The 
steady-state  hover  condition  results  in  the  following 
simplifications  of  the  nonlinear  system  equations. 

1.  The  product  of  two  small  numbers  is  an  extremely  small 
number,  thus  terms  involving  the  products  of 
translational  or  angular  velocities  are  equal  to  zero 
(e.g.  rq,  pq,  (pw-ru)  are  set  equal  to  zero) . 

2.  The  aerodynamic  forces  (F„,  F,y,  Fu)  and  moments  (L,,  M,, 
NJ  are  very  small  for  a'  equal  to  90°,  and  can  be 
neglected  in  the  hover  flight  condition.  (Note: 
Because  the  vanes  lie  within  the  downwash  of  the 
propeller,  the  vane  effectiveness  coefficients  can  not 
be  neglected.) 

3.  The  sine  of  a  state  is  equal  to  the  state  and  the 
cosine  of  a  state  is  equal  to  one.  This  is  the  small 
angle  approximation  for  angles  less  than  15  degrees. 

4.  is  a  small  number  equal  to  0.0729.  For  hover  or 
very  small  translational  velocities,  5rpia  is  a  very 
small  number.  Therefore,  the  product  of  and  6rpm 
is  neglected. 

5.  The  propeller  angular  velocity,  a>p,  is  considered 
constant.  Thus,  the  propeller  angular  acceleration, 
wp,  is  equal  to  zero. 


Table  2 . 3  lists  the  Archytas  system  equations  when  the  above 
simplifications  are  applied  to  the  nonlinear  system 
equations.  Note  that  much  of  the  coupling  between  states 
and  all  of  the  nonlinear  products  of  states  have  been 
eliminated. 
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Table  2.3  ARCHYTAS  LINEARIZED  HOVER  EQUATIONS 


The  force,  F-p^,  in  Table  2.3  is  a  function  of  the 
engine  rpm  and  is  computed  by  equation  (2.46).  The  moment 
terms,  Lt,  M,  and  Nt,  are  functions  of  the  displacement  of  the 
aileron,  elevator  and  rudder.  They  are  computed  by  equation 
(2.48).  The  moments  of  inertia,  Ix,  Iy,  Iz  and  I„,  were 
determined  by  the  AROD  project  engineers  and  are  listed  in 
Appendix  A.  The  angular  velocity  of  the  propeller,  wp,  in 
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the  hover  mode  equals  712.0943  rad/sec  (6800  rev/min) .  The 
weight  of  the  AROD  prototype  is  equal  to  76.5  lbs.  Using 
these  constant  parameters,  the  forces  and  moments  for  the 
hover  steady  state  reference  condition  were  computed  and  are 
listed  in  Table  2.4.  Substituting  the  results  of  Table  2 . 4 
into  the  equations  of  Table  2 . 3  yields  the  linearized  hover 
equations,  which  are  summarized  in  Table  2.5. 

D.  SUMMARY 

The  nonlinear  differential  equations  of  motion  of  a 
rigid  body  were  developed  from  Newton's  second  law.  The 
equations  were  then  modified  for  the  specific  case  of 
Archytas.  This  modification  included  the  effect  of  the 
spinning  three-bladed  rotor.  The  equations  were  linearized 
about  a  steady-state  hover  condition.  Next,  the  applied 
forces  and  moments,  and  the  servo  and  engine  equations  were 
defined.  The  forces  and  moments  specific  for  the  hover 
reference  condition  were  computed  and  applied  to  the 
Archytas  linear  model. 
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Table  2.4  FORCES r  MOMENTS  AND  CONSTANTS 


Moments  Forces 

Lt=Laeffha  FThtx=  THOVER  *  XRPM  5„, 

£,,,=  -150,379 . 57  5  ^  FTAr=76 . 5+0 . 23  876rpB 

factor^  e 

Mt= (-112,716 .87)  (1)6, 

~Nte££Y£ actor1 &  r 

^(.=  (-128,774.80)  (l)6r 

Moments  of  Inertia 

Jx=7  06  3 . 39  ,  Jy=7768.22,  Iz=7729.58,  Jrx=69.552  {Ib^  in2) 

Angular  velocity 

o>  =712.0943-^; 

*  sec 

Weight /gravity /mass 

weights  6. 5  lbs.}  gravity*  3  2. 174-^-;  /nass= 

sec2  gravi ty 

Engine  Lag  Time  Constant  and  Scale  Factor 
Wj?=2.0  837.758  rad/ sec/ rad 
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Table  2.5  LINEARIZED  HOVER  STATE  EQUATIONS 


SI&XE 


a=  “150,379 .57  *  _  _ 

P - 7~ 6~3.39  ^-21.295, 


g= 


r= 


7768.22 

1 

7729.58 


14 .514,-6 .75r 


[-112,716 .875^-52,440. 82r]  = 

[-12  8,774.805-  +  52,440.82g] =-16 . 685r+6 .78 g 


.  (76.5+0. 2387  6„_) 

U= - 2^8 - ^“32.174=0.100295^ 

V-=  32.174? 

V=-32 . 1740 


6=g 

4f=p 

?=r 

*as-«i  *a~H2&a+H2Ua 

fir=’^i4r-^25r^2Ur 

S  t  ^1^  t-^2®  t+^2^t 
^  rpm  ~  “  01 rpm + W  t 
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III.  OPTIMAL  CONTROL  THEORY 

The  purpose  of  this  chapter  is  to  develop  a  procedure 
for  applying  optimal  control  theory  to  the  solution  of 
tracking  problems.  First,  several  reasons  for  desiring  an 
optimal  solution  are  presented.  Next,  the  state  space 
representation  (both  continuous  and  discrete)  of  a  system  is 
given.  Finally,  the  application  of  optimal  control  to  the 
solution  of  tracking  problems  is  illustrated. 

A.  WHY  OPTIMUM  CONTROL  FOR  ARCHYTAS? 

The  first  reason  for  seeking  an  optimum  controller  is 
that  feedback  gains  can  be  computed  for  a  much  broader  range 
of  control  problems.  Specifically,  optimal  control  provides 
solutions  for  high  order,  multiple-input,  multiple-output 
(MIMO)  systems.  Such  systems  are  often  intractable  with 
classical  methods.  The  pitch  and  yaw  angle  controller  for 
Archytas  is  an  eight  order  MIMO  system.  Thus,  optimal 
control  is  the  preferred  method. 

Additionally,  optimal  control  lends  itself  nicely  to  a 
discrete  time  solution  of  the  control  problem.  Archytas 
will  employ  an  on  board  digital  computer  to  perform  inflight 
stability  and  control.  While  a  continuous  time  controller 
can  be  easily  discretized  in  many  cases,  design  of  a  sampled 
data  controller  will  simplify  the  procedure. [Ref .  3:  p.  58] 
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Finally,  optimal  control  is  known  to  provide  robust  and 
insensitive  solutions  to  the  feedback  control  problem. 
Assuming  that  an  appropriate  performance  measure  is  chosen 
to  determine  the  optimal  feedback  gain  matrix,  K,  the 
solution  can  be  expected  to  have  a  fair  degree  of  tolerance 
to  plant  model  inaccuracies.  Clearly,  robustness  is  not 
only  a  desired  property  of  the  controller,  it  is  an  absolute 
necessity  if  the  controller  is  to  be  applicable  to  both  the 
linear  and  nonlinear  models  of  Archytas. [Ref .  3:  pp.  58-59] 

B.  STATE  SPACE  REPRESENTATION 

The  state  of  a  system  may  be  defined  as  the  minimum 
amount  of  information  required  such  that  (given  the  input  to 
the  system)  the  response  of  the  system  is  completely 
determined  for  all  future  time.  For  dynamic  systems,  the 
response  of  the  system  is  defined  by  the  differential 
equations  that  model  the  system,  the  initial  conditions  and 
the  forcing  function.  The  number  of  state  variables  or 
states  is  equal  to  the  total  order  of  the  systems 
differential  equations.  In  order  to  provide  a  systematic 
mathematical  approach  to  analysis  of  the  system,  it  is 
convenient  to  describe  the  system  by  a  set  of  first-order 
differential  equations.  This  set  of  equations  is  called  the 
state  equations  and  constitute  the  basis  for  the  state  space 
representation.  [Ref.  8:  pp.  2C6-207] 
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l.  Continuous  Tina  Systems 

The  state  space  representation  of  a  general  n*  order 
continuous -time,  time-invariant  system  is  described  by  the 
following  matrix  state  equations: 

Kit)  =4  s(e)  u(t)  ;  (3.1) 

£(t)=£2C(e)  ;  (3.2) 

where  the  definitions  in  Table  3.1  apply  to  a  system  with  £ 
control  inputs  and  m  measurable  outputs.  The  equations 

TABLE  3.1 


STATE  SPACE  DEFINITIONS  FOR  CONTINUOUS -TIME  SYSTEMS 


Term 

Dimension 

Definition 

Kit) 

in  x  1) 

State  vector 

Xit) 

{m  x  l) 

Output  vector 

4 

(n  x  n) 

Plant  matrix 

M 

in  x  £) 

Control  distribution  matrix 

£ 

im  x  n) 

Output  Distribution  matrix 

listed  in  Table  2.5  are  linear  and  t ime- invariant . 

2.  Discrete  Time  Systems 

As  was  noted  earlier,  there  are  many  benefits  for 
seeking  a  discrete  time  solution  for  the  Archytas  control 
problem.  Therefore,  the  automatic  control  systems  designed 
will  focus  on  the  application  of  optimal  control  theory  to 
discrete  time  systems. 
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Similar  to  the  continuous-time  system,  the  state 
space  representation  of  a  general  n111  order  discrete-time 
system  is  described  by  the  following  matrix  state  equations: 

2f(/c+l)  =<£  x(k)  +£  Ji(&)  ;  (3.3) 

£(k)=£x(k)  ;  (3.4) 

JL  and  E  are  defined  as: 


£=e*c  ;  (3.5) 

£=fQTe*cdt  B  ;  (3.6) 

where  T  is  the  sampling  period  and  k  is  an  integer  time 
index.  Reference  10  provides  a  detailed  development  of  the 
relationship  between  continuous-time  and  discrete-time 
systems.  The  definition  in  Table  3.2  apply  to  a  system  with 
l  control  inputs  and  m  measurable  outputs. 

TABLE  3.2 


STATE  SPACE  DEFINITIONS  FOR  DISCRETE-TIME  SYSTEMS 


Term 

Bimsaaion 

Detinition 

X(k) 

(n  X  l) 

State  vector 

xik) 

{m  x  l) 

Output  vector  (0  x  m  s  n) 

a 

(n  x  n ) 

Plant  matrix 

£ 

(n  x  H) 

Control  distribution  matrix 

Q 

( m  x  n) 

Output  Distribution  matrix 
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c.  LINEAR  QUADRATIC  REGULATOR  PROCEDURE 

The  theory  behind  the  linear  quadratic  regulator  is  well 
developed.  Procedures  exist  which  make  the  design  of 
controllers  using  linear  quadratic  regulator  theory  easily 
obtainable.  The  purpose  of  this  section  is  to  illustrate 
the  application  of  the  linear  quadratic  regulator  to  the 
tracking  problems  associated  with  Archytas.  Reference  11 
provides  a  detailed  development  of  the  linear  quadratic 
regulator  theory. 

1.  Quadratic  Cost  Function 

The  discrete  LQR  synthesis  problem  is  that  of 
determining  the  control  that  minimizes  the  performance 
measure: 


J^xT{k)Q  x(k)  *  uT(k)  R  u(k)  ;  (3.7) 

k»  o 


where 

J  =  Scalar  cost  of  operating  the  system; 

x(k)  =  State  vector  at  discrete  times; 

u(k)  =  Control  vector  at  discrete  times; 

£>  =  Symmetric  state  trajectory  weighting 

matrix; 

R  =  Symmetric  control  weighting  matrix; 

T  =  Matrix  transpose  operator. 
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The  solution  to  this  optimization  problem  is  the  linear 
controller: 

u(k)  =- (M~1£T  tf)x(Jc)  =~KX(k)  ;  (3.8) 

where  M  satisfies  the  algebraic  Riccati  equation  (ARE) : 

&  =  +  &TM  ~  MM  BrlgTU  +  £  •  (3.9) 

[Ref.  11:  pp.  345-346 ]  Many  software  packages,  including 
the  program  MATLAB  used  in  this  thesis,  contain  subroutines 
that  calculate  the  value  of  K  for  a  given  dynamic  system  and 
performance  measure. 

2.  Performance  weighting  Matrices 

The  optimal  control  is  fully  specified  by  the  system 
model  and  the  weighting  matrices  2  an<*  £•  2  penalizes 

deviation  of  the  state  vector  £  from  the  origin  and  R 
penalizes  the  use  of  too  much  control  effort.  Trial  and 
error  was  used  in  selecting  values  for  these  performance 
weighting  matrices.  The  relationship  between  the  weighting 
matrices  2  and  R  and  the  dynamic  behavior  of  the  closed-loop 
system  depends  of  course  on  the  matrices  A  and  IJ  and  is 
quite  complex.  The  approach  taken  in  the  design  of  the 
controllers  for  Archytas  was  to  solve  for  the  gain  matrices 
2  that  result  from  a  range  of  weighting  matrices  £  and  R, 
and  then  simulate  the  corresponding  closed-loop  response. 

The  gain  matrix  K  that  produced  the  response  closest  to 
the  desired  design  objectives  was  chosen.  [Ref.  11:  p.  341] 
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3.  Optimal  Tracking  Systems 

The  regulator  and  tracking  problems  appear  very 
similar.  The  tracking  problem  attempts  to  drive  the  states 
of  the  system  to  a  desired  level;  whereas,  the  regulator 
attempts  to  drive  all  of  the  states  to  zero.  Their 
differences  present  conceptual  difficulties  and  sensitivity 
problems  when  viewed  in  a  practical  context  [Ref.  12:  pp, 
643-647],  Burl  [Ref.  13]  presents  a  comprehensive 
development  of  the  subtleties  encountered  when  applying  LQR 
synthesis  to  the  tracking  problem.  This  development  is 
generalized  below  to  demonstrate  these  subtleties  and  to 
indicate  a  procedure  for  selecting  the  proper  form  of  the 
performance  measure.  The  following  development  is  applied 
to  a  first  order  system,  but  the  results  are  applicable  to 
systems  of  arbitrary  order. 

Given  the  scalar  plant: 

y{  t)  =~Ay(  t)  +Bu(  t)  .  (3.10) 

The  purpose  of  the  control  system  is  to  drive  the  output 
y(t)  to  the  constant  reference  input  r.  This  results  in  the 
error  equation: 

e(t)=z-y(t)  ;  (3.11) 

where  the  desire  is  to  drive  e(t)  to  zero.  The  application 
of  linear  regulator  theory  requires  that  this  error  be  a 
linear  combination  of  the  states  of  the  system.  This  can  be 
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readily  accomplished  by  generating  a  new  state  equation  for 
e(t)  : 


e{  t)  =-y(  t)  =Ay  ( t)  ~Bu  ( t)  ; 

adding  and  subtracting  Ar  yields: 

e(t)  =Ay{t)  +Ar-Ar-Bu(t)  ; 

6(t)  =A[-r+y(t)]  -Bu(t)  +Ar  ; 
e(t)  =~Ae(t)  -Bu(t)  +Ai  . 

Since  the  error  equation,  Equation  (3.19),  is  linear  and 
time  invariant,  the  performance  measure: 

00 

J=f  [Qe2(t) +u2(t)}dt  ;  (3.14) 

o 

should  provide  an  optimal  solution  for  the  system  of 
Equation  (3.21).  However,  as  shown  by  Burl  [Ref.  13]  this 
optimal  control  problem  will  not  yield  a  solution.  To 
demonstrate  this  fact,  note  that  the  existence  of  this 
integral  requires  that  both  u(t)  and  e(t)  approach  zero  as  t 
tends  to  infinity.  If  e(t)  approaches  zero,  then  e ( t)  must 
also  approach  zero  which  from  Equation  (3.21)  yields: 

0  =  -A ( 0 )  -Bu{t)  +Ar  -  U{t)  *4r  *  0  (3.15) 

& 

This  nonzero  value  for  u(t)  implies  that  the  performance 
measure  is  infinite.  A  solution  to  this  problem  would  be  to 
formulate  the  performance  measure  as 


(3.12) 


(3.13) 
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(3 . 16) 


J=j [QeHt)  +{u(t)  ~^r)2]dc  . 


This  will  result  in  the  existence  of  a  theoretically  optimal 
control  provided  that  the  model  of  the  system  (the 
parameters  A  and  B)  is  exactly  known.  When  errors  exist  in 
the  model,  a  nonzero  steady  state  error  will  exist.  The 
resulting  control  system  will  be  mathematically  correct,  but 
it  will  be  unacceptable  for  many  applications  due  to  its 
sensitivity  to  changes  in  the  plant  parameters.  The  control 
system  will  be  super-tuned  (it  will  not  be  robust) . [Ref .  12: 
p.  645] 

This  undesirable  result  can  be  overcome  by  letting  the 
controller  find  the  appropriate  steady  state  value  of  the 
control.  This  is  achieved  by  application  of  the  performance 
measure : 


J=f[Qe2(t)  +ii2{t)]dt  .  (3.17) 

0 

The  control  that  minimizes  Equation  (3.25)  given  the  system 
of  Equation  (3.21)  is  found  by  first  differentiating 
Equation  (3.21),  yielding: 

e(t)  =-A6(t)  -Bill  t)  .  (3.18) 

Noting  that  the  error  is  the  integral  of  £(t)  results  in 
the  state  space  system: 
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(3.19) 


£(  t) 

Jo  1 ' 

e{  t) 

0  ' 

e{  c) 

JO  -A 

,<Mt)J 

J-3 

This  system  and  the  performance  measure  of  Equation  (3.25) 
form  a  linear  regulator  problem  whose  solution  is  state 
feedback: 


ti(t)  =-  [k. 


k2] 


e(t) 

£(t). 


(3.20) 


The  gains  k,  and  k2  are  computed  by  application  of  linear 
regulator  theory.  The  actual  control  that  is  applied  to  the 
system  is  found  by  integrating  Equation  (3.28): 


t 

u(  t)  =~k1fe(x)  dt-k2e(  fc)  (3.21) 

o 

The  resulting  controller  incorporates  proportional  plus 
integral  feedback  of  the  error.  This  fact  is  reasonable 
since  the  system  to  be  controlled,  Equation  (3,21),  is  of 
type  o  with  a  steady  state  disturbance  input. 

To  summarize,  when  LQR  synthesis  is  applied  to  the 
tracking  problem,  the  proper  choice  of  state  variables  will 
help  to  eliminate  sensitivity  problems  and  ensure  system 
robustness.  The  state  variables  should  not  contain  any 
input  dependent  terms,  and  they  should  asymptotically 
approach  zero.  The  system  tracking  error  and  its  time 
derivative  are  natural  state  variable  for  the  LQR  synthesis 
procedure.  In  essence,  the  key  to  a  successful  design 
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depends  upon  the  proper  formulation  of  the  performance 
measure.  The  performance  measure  of  Equation  (3.25) 
represents  one  possible  formulation  that  yields  reliable 
results.  This  performance  measure  is  used  in  the  design  of 
the  Archytas  control  systems.  These  ideas  are  extended  to  a 
general  multiple  input,  multiple  output  (MIMO)  system  by 
Burl  [Ref .  13 ] . 

D.  SUMMARY 

In  this  chapter,  the  application  of  the  linear  quadratic 
regulator  to  the  tracking  problem  has  been  addressed.  A 
procedure  for  formulating  a  performance  measure  applicable 
to  tracking  problems  was  developed.  Additionally,  the  state 
space  representations  for  both  continuous-time  and 
discrete-time  systems  were  presented.  In  the  next  chapter, 
the  results  of  this  chapter  are  applied  to  generate  control 
systems  for  Archytas. 
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IV.  CONTROL  SYSTEM  DESIGN  FOR  ARCBYTAS 

The  purpose  of  this  chapter  is  to  use  optimal  control 
theory  to  design  an  automatic  flight  control  system  for 
Archytas  during  the  hover  mode  of  flight.  Because  of  the 
coupled  nature  of  the  Archytas  control  problem,  a  linear 
quadratic  regulator  (LQR)  synthesis  technique  is  pursued. 

A.  ARCBYTAS  CONTROL  SUBSYSTEMS 

The  automatic  control  system  is  logically  separated  into 
three  subsystems  according  to  the  linearized  equations  of 
Table  2.5.  The  three  control  subsystems  are: 

1.  Roll  rate  controller; 

2.  Altitude  rate  controller; 

3.  Pitch  angle  and  yaw  angle  controller. 

Because  each  of  these  control  subsystems  is 
independently,  any  cross-coupling  which  may 
the  subsystems  is  not  taken  into  account. 

B.  ARCHYTAS  ROLL  RATE  CONTROLLER 
1.  The  Roll  system 

The  roll  rate  controller  for  Archytas  will  serve 


designed 
occur  between 


three  primary  functions: 


1.  Allow  the  operator  to  command  a  desired  rotation 
velocity  about  the  vehicle's  longitudinal,  or  x,  axis. 
This  rotation  velocity  will  permit  the  operator  to 
position  Archytas  by  terminating  the  roiling  motion 
when  a  desired  angle  is  achieved.  This  capability  is 
critical  during  the  landing  phase  in  order  to  position 
Archytas  correctly  with  respect  to  the  wind. 

2.  Eliminate  the  rotation  velocity  imparted  to  Archytas 
from  the  effect  of  cross-winds. 

3.  Eliminate  the  rotation  velocity  imparted  to  Archytas 
from  the  reactive  torques  applied  to  the  roll  axis  as 
the  engine  speed  is  varied. 


The  simplified  equation  of  motion  which  describes 
Archytas'  roll  rate  subsystem  is  given  as: 

p=-21 , 296a  .  (4.1) 

The  aileron  servo  dynamics  are  modeled  in  Appendix  B  as  a 
second  order  dynamical  system  with  a  natural  frequency,  ojn, 

of  52.4  rad/sec  (8.34  Hz)  and  a  damping  ratio,  £,  of  0.707. 
The  corresponding  differential  equation  is: 


5a=-74.l6a-2745.85a+2745.8ua  .  (4.2) 

From  Equations  (4.1)  and  (4.2),  the  state  equation 
can  be  written  in  the  matrix  form  &=&2£+||u: 


p 

0  -21.29  0 

P 

■  0  ■ 

5a 

K 

0  0  1 

0  -2745.8  -71.1 

+ 

0 

2745 . 8. 

(4.3) 


The  roll 


system  tracking  error  is  defined  as: 
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Ep~pc~p  ; 


(4.4) 


where  pc  is  the  input  command  and  p  is  the  measured  roll 
rate.  From  Equation  (4.4),  the  differential  equation  for 
the  tracking  error  is: 

4='P  .  (4.5) 

If  Equation  (4.5)  is  combined  with  the  time  derivatives  of 
Equations  (4.1)  and  (4.2)  a  new  state  equation  can  be  formed 
that  is  appropriate  for  tracking  system  design  (as  discussed 
in  Chapter  III) .  This  state  equation  is: 


*p 

0 

-1 

0 

0  ' 

'EP 

0 

p 

0 

0 

-2a . 29 

0 

0 

K 

0 

0 

0 

1 

■f 

0 

s. 

0 

0 

-2745.8 

-74.1 

K 

27  45.8. 

2.  Roll  Rate  controller  Design 

a.  Sampling  Frequency  Selection 

The  application  of  optimal  control  theory  to 
discrete-time  systems  requires  that  an  appropriate  sampling 
frequency  be  determined.  The  sampling  frequency,  fs,  is 
simply  the  inverse  of  the  sampling  period,  T.  A  general 
rule  of  thumb  in  control  systems  is  to  sample  a  system  such 
that 
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(4.7) 


fs*1 04ax  / 

where  is  the  maximum  bandwidth  of  the  system.  From 
Appendix  B,  the  natural  frequency  of  the  servos  is  52.4 
rad/sec  (8.34  Hz).  Bassett  [Ref.  3:  p.  76]  demonstrates  how 
the  Archytas  MIMO  system  can  be  considered  as  a  set  of 
several  SISO  systems.  Each  system  has  a  different  bandwidth 
for  the  purpose  of  determining  a  sampling  frequency.  The 
subsequent  highest  natural  frequency  is  equal  to  4.64 
rad/sec  (0.74  Hz).  Because  the  natural  frequency  of  the 
servos  is  a  factor  of  ten  greater  than  the  highest  natural 
frequency  of  4.64  rad/sec,  the  servo  natural  frequency  will 
be  used  to  determine  the  required  sampling  period,  T.  From 
Equation  (4.7),  f =10(52.4  rad/sec) =524 . 0  rad/sec  (83.4  Hz), 
using  83.4  Hz  as  the  sampling  rate  yields 

t=T  =  '8TT=0'012  sec  *  (4,8> 

For  the  controller  designs  of  this  thesis,  a  sampling  period 
of  0.01  seconds  will  be  used. 

b.  Discretizing  the  Roll  Rate  System 

MATLAB  provides  tools  for  computing  the  discrete¬ 
time  matrix  state  equation.  With  the  sampling  period  T 
equal  to  0.01  seconds,  the  discrete-time  state  space  can  be 
written  in  the  matrix  form  x(k+l)  =£  x(k)  +£  u(k)  : 
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Bp(k+lV 

'1 

-0.01  0.001 

p{k+i) 

0 

1  -0.2048 

&a(k+i) 

0 

0  0.8935 

S4(k+ l) 

0 

0  -18.5258 

0  • 

■Ep(k)' 

0 

-0 . 0008 

p(k) 

-0 . 0081 

0.0067 

&s<k> 

+ 

0.1065 

0 . 3936  . 

18 . 5258. 

(4.9) 


c.  Gain  Determination 

The  optimal  control  is  determined  from  the  state 
equation  and  the  performance  measure: 


Bp{kV 

[Ep(k)  p(k)  6 a(k)  Sa(k)]£ 

p{k) 

*a(k) 

5.1*1, 

+  uT(k)g  u{k) 

Jt»0 


The  optimal  control  u(k)  is: 


(4.10) 


**(*)  *  "  [*i  K  k3  k4] 


'Bp(k) 

p(k) 


•-&£(k)  (4.11) 


where  K  is  the  steady-state  gain  matrix.  The  actual  control 
u(k)  is  then  obtained  by  integrating  u{k)  to  obtain: 


u(k)  =£  u(w)  =  I*!  k2  k 3  k4] 


Ep(m)] 
P(m) 
6a(m) 
&a{n>) 

ep^  ~k2P(™)  ~K36a(m)  -k4ta  (m) 


.11*0 


.11*0 


k 

L 

w-0 


(4.12) 
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The  LQR  weighting  matrices,  2  and  E'  are  chosen  to  satisfy 
the  following  design  criteria: 

1.  The  overshoot  to  a  step  input  should  be  less  than  five 
percent. 

2.  The  five  percent  settling  time,  t5„,  is  less  than  or 
equal  to  1  second. 


Using  an  iterative  process,  the  weighting  matrices  that  meet 
the  above  design  criteria  were  found  to  be: 


.3  0  0  0 
0  0  0  0 
0  0  0  0 
0  0  0  0. 


S=[l]  • 


(4.13) 


The  resulting  steady-state  feedback  gain  matrix,  K,  is: 


Kr  IK  ^  *3  *41  *  [0.5347  -0.2385  0.1329  0.0017]  .  (4.14) 

Figure  4.1  shows  the  roll  rate  controller  block  diagram. 
d.  Simulation  Results 

Figure  4.2(a)  shows  the  response  of  the  closed 
loop  system  to  a  step  input  of  six  degrees /second.  The 
design  criteria  of  an  overshoot  less  than  five  percent  and 
the  five  percent  settling  time,  t5„,  less  than  one  second 
are  achieved.  Figure  4.2(b)  shows  the  aileron  vane  angle  as 
a  result  of  the  step  input.  Weir  [Ref.  16]  demonstrates 
that  the  control  vanes  stall  when  displaced  by  an  angle  of 
plus  or  minus  30  degrees  from  the  air  flow  zero  reference. 
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Therefore,  it  is  important  to  ensure  that  the  design  not 
require  vane  angles  greater  than  30  degrees. 

Because  the  roll  rate  controller  was  designed 
using  the  linear  model,  it  is  necessary  to  measure  the 
relative  robustness  of  the  design  to  ensure  that  it  will  be 
applicable  to  the  nonlinear  model.  The  phase  and  gain 
margins  are  such  measures.  From  Figure  4.3,  the  gain  margin 
is  equal  to  32.63  dB  and  the  phase  margin  is  equal  to  63.99 
degrees.  A  general  rule  of  thumb  is  to  require  a  gain 
margin  greater  than  six  dB  and  a  phase  margin  greater  than 
30  degrees  to  ensure  satisfactory  performance.  Clearly, 
these  benchmark  values  are  exceeded  and  the  above  design 
should  perform  well  when  applied  to  the  nonlinear  mod  1. 

The  steady-state  gains  of  Equation  (4.14)  are  an 
optimal  solution  for  the  performance  measure  of  Equation 
(4.10).  However,  this  optimal  solution  requires  that  all  of 
the  states  be  measured  or  estimated.  Although  measurement 
of  the  vane  angle,  5a,  is  possible,  a  computational  observer 

would  have  to  be  included  to  provide  an  estimate  of  the  vane 
velocity,  S3.  Inclusion  of  an  observer  would  add  an 

additional  task  to  the  onboard  digital  computer,  and 
increase  the  complexity  of  the  controller  design.  An 
alternative  to  the  observer  would  be  simply  to  set  the  vane 
velocity  gain  to  zero.  The  system  would  have  to  be 
simulated  with  this  gain  equal  to  zero  and  the  phase  and 
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gain  margins  computed  to  determine  how  the  system  response 
would  be  affected. 

e.  Reduced  Order  Model 

Because  the  servo  dynamics  are  significantly 
faster  than  all  other  system  dynamics,  a  second  alternative 
is  to  ignore  the  servo  dynamics  completely  and  form  a 
reduced  order  state  space  representation  of  the  system. 


The  discrete-time  state  equation  is  given  as 


(4.15) 


•Ep{Jc+l)l  [l.o  -O.Oll  [V*>|  +  [  0.0011  ]  6  (4. 16) 

.£<*+1)  [o  l.oJ[£(lc)J  1-0.2129]  a 


Now,  the  steady-state  gain  matrix,  K,  is  determined  from 
Equation  (4.16)  and  the  performance  measure: 


=  g  jiV*)  pm  jj2  ^  *  a^ms  o(«j  .  <4-17’ 


The  optimal  control  u(k)  is: 


\Ea{k) 

am  ■  -(*z*2lLfc) 


=  -K &{k) 


(4.18) 


and  u(k)  is 
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(4.19) 


u[k)  =  £  awn)  =  -]P  [kz  kz] 

)  .-!*0 
k 

■k.  ^  E?(m)  -  kzp{m) 


.71*0 


K„(/n)l 

i 

LP(^)  J 


The  weighting  matrices,  £>  and  R,  are  defined  as 


0.3  0 
0  0 


£=[13 


(4.20) 


The  resulting  steady-state  feedback  gain  matrix,  K,  is 

£=(0.4376  -0.2027]  .  (4.21) 

Figure  4.4(a)  shows  the  response  of  the  system  to 
a  step  input  of  six  degrees /second.  The  design  criteria  for 
overshoot  and  settling  time  are  achieved.  The  step  response 
for  this  reduced  order  model  is  almost  identical  to  the 
response  for  the  full  order  model  given  in  Figure  4.2(a). 
From  Figure  4.4(b),  the  maximum  magnitude  of  the  aileron 
vane  angle  displacement  is  0.47  degrees.  This  vane 
displacement  angle  meets  the  requirement  of  less  than  or 
equal  to  30  degrees. 

The  necessity  of  determining  the  relative 
robustness  of  the  system  is  now  magnified  due  to  the  reduced 
order  controller  design.  From  Figure  4.5,  the  gain  and 
phase  margins  are  equal  to  20.59  dB  and  57.66  degrees 
respectively.  Note  that  these  margins  were  computed  for  the 
controller  with  the  full  order  model.  The  reduction  in  gain 
and  phase  margins  from  the  full  order  model  is  an  expected 
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Figure  4. 4  deduced  order  Controller 
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Figure  4.5  Reduced  Order  Controiler 
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consequence  of  the  reduced  order  controller.  However,  the 
requirement  for  a  minimum  of  six  dB  gain  margin  and  20 
degrees  of  phase  margin  is  satisfied.  Therefore,  it  is 
expected  that  the  reduced  order  controller  design  will  be 
applicable  to  the  nonlinear  model.  Figure  4.6  is  a 
graphical  realization  of  the  reduced  order  controller 
applied  to  the  roll  rate  system. 
f .  Summary 

An  optimal  tracker  was  designed  for  the  roll  rate 
control  by  applying  the  procedures  of  chapter  III.  This 
control  system  yielded  very  satisfactory  system  performance 
and  robustness.  Recognizing  the  dynamics  of  the  servos 
could  possibly  be  ignored,  a  second  reduced  order  design  was 
developed.  The  resulting  controller  proved  to  be  robust  and 
almost  identical  in  performance  to  the  full  order  controller 
design. 

Because  of  the  favorable  results  obtained  with 
the  reduced  order  roll  rate  controller  design,  the  altitude 
rate  and  the  pitch  and  yaw  angle  controllers  will  be 
designed  using  this  technique.  The  resulting  altitude  rate 
and  pitch  and  yaw  controllers  will  be  evaluated  analogously 
to  the  roll  rate  controller  to  ensure  they  meet  or  exceed 
the  design  criteria. 
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c.  ARCHYTAS  ALTITUDE  RATE  CONTROLLER 
1.  The  Altitude  System 

The  primary  function  of  the  altitude  rate  controller 
for  Archytas  is  to  allow  the  operator  to  command  a  desired 
translational  velocity  along  the  vehicled  longitudinal,  or 
x  axis.  This  translational  velocity  will  permit  the 
operator  to  position  Archytas  at  a  given  altitude  by 
terminating  the  velocity  when  a  desired  altitude  is 
achieved.  This  capability  is  critical  during  the  landing 
phase  in  order  zo  control  the  altitude  as  Archytas  lands  in 
a  vertical  position. 

The  simplified  equation  of  motion  which  describes 
Archytas'  altitude  rate  subsystem  is  given  as: 

U  =  h  =  0.10029  6rp„  (4.22) 

Note  the  change  of  notation  from  u  to  h ,  which  is  more 
appropriate.  The  throttle  servo  model  is  identical  to  the 
aileron  servo  model  and  is  given  as: 

6t=-74  .  l6t-2745  .  S6t+2745  .  Qu,  .  (4.23) 

The  Archytas  engine  is  modeled  as  a  first  order  lag: 


-^.,+  167  5 


(4.24) 


The  altitude  rate  tracking  error  is  defined  as: 
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E^he-h  ; 


(4.25) 


where  Ac  is  the  input  command  and  h  is  rhe  measured 

altitude  rate.  From  Equation  (4.25),  the  differential 
equation  for  the  tracking  error  is: 

E^-h  .  14*26 

Combining  Equation  (4.26)  with  the  time  derivatives  of 
Equations  (4.22),  (4.23)  and  (4.24),  a  state  equation  that 

is  appropriate  for  tracking  system  design  is: 


0  -1 
0  0 
0  0 
0  0 
0  0 


0 

-0.10029 

-2 

0 

0 


0 

0 

1675.5 

0 

-2745.8 


(4.27) 


In  a  manner  identical  to  that  used  for  the  roll  rate 
controller  design,  the  reduced  order  altitude  rate  model  is 
formed.  The  controller  is  designed  using  the  reduced  order 
model  and  its  performance  and  robustness  are  evaluated  with 
respect  to  the  full  order  model.  Ignoring  the  servo 
dynamics,  the  reduced  order  state  space  model  is: 


0  *6  0 

0.10029  h  +  0  ut  . 

-2  6  1675.5 


(4.28) 
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2.  Altitude  Rate  controller  Design 

a.  Discretizing  the  Altitude  Rate  System 
With  the  sampling  period,  T,  equal  to  0.01 

seconds,  the  discrete-time  state  equation  is: 

Ef,(k+1)  1-0.01  o  E&(k)  j-  0  I 

h[k+i)  =o  1  0.001  h{k)  +|  0  j  u;  .  (4*29) 

iIpmik*u\  [o  0  0.9802J  [l675 . SJ 

b.  Gain  Determination 

The  LQR  weighting  matrices,  Q  and  S'  for  the 
altitude  rate  control  loop  are  chosen  to  satisfy  the 
following  design  criteria: 

!•  The  overshoot  to  a  step  input  should  be  less  than  five 
percent. 

2.  The  five  percent  settling  time,  t5%,  is  less  than  or 
equal  to  2  seconds. 

The  resulting  weighting  matrices  are: 

0.1  0  o' 

Qr  0  0.01  0  ;  £»[l]  .  (4.30) 

0  oo. 

The  steady-state  feedback  gain  matrix,  X,  is: 

&-  [-0.3060  0.2003  0.0038]  .  (4.31) 

Figure  4.7  shows  the  altitude  rate  controller  block  diagram. 
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Figure  4.7  Reduced  Order  Altitude 
Rate  Block  Diagram 


c.  Simulation  Results 

Figure  4.2(a)  shews  the  response  of  the  altitude 
rate  closed  loop  system  to  a  step  input  of  ten  feet/ second. 
The  design  criteria  of  an  overshoot;  less  than  five  percent 
and  the  five  percent  settling  time,  t5<t ,  less  than  two 
seconds  is  achieved.  Figure  4.3(b)  shows  the  control 
applied  to  the  system. 

The  gain  and  phase  margins  were  determined  as  a 
measure  of  the  system  robustness.  From  Figure  4.9,  the  gain 
margin  is  equal  to  18.03  aB  ana  the  phase  margin  is  equal  to 
58.78  degrees.  These  values  of  gain  and  phase  margin 
indicate  that  the  altitude  rate  control  system  design  based 
on  the  reduced  order  model  should  provide  good  results  when 
applied  to  the  nonlinear  system  model. 

D.  ARCHYTAS  PITCH  AND  YAW  ANGLE  CONTROLLER 
1.  The  Pitch  and  Yaw  Angle  System 

The  purpose  of  the  pitch  and  yaw  angle  controller 
for  Archytas  is  to  maintain  commanded  orientation  around  the 
pitch  and  yaw  axes.  This  requirement  is  necessary  to  allow 
the  operator  to  position  Archytas  during  landing  ov  pitching 
or  yawing  the  vehicle  slightly  to  induce  a  translationai 
velocity  along  the  body  fixed  y  or  z  axis.  This  control 
problem  is  complicated  by  the  gyroscopic  coupling  between 
the  axes  caused  by  the  propeller. 
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Reduced  Order  Discrete  Open  Loop  Bode  Plot 


Figure  4.9  Reduced  Order  Altitude  Rate 
Controller  Gain  and  Phase  Margins 


gp  uibo 


The  simplified  equations  of  motion  that  describe  the 
pitch  and  yaw  angle  subsystem  are  given  as: 

£=-11 . 196e-i . 52r;  (4.32) 

f=-12.285r+1.78g  .  (4.33) 

The  elevator  and  rudder  servos  are  modeled  in  a  manner 
identical  to  the  aileron  servo  and  the  corresponding 
differential  equations  are: 


and 


6e=-74.l6e-2745.86e+2745.Sue  ;  (4.34) 


5r=-74.l6r-2745.88r+2745.8ur  .  (4.35) 

The  pitch  and  yaw  angle  tracking  errors  are  defined 

as: 

i?e=0c-e  ;  (4.36) 

and 

;  (4.37) 

where  0C  and  Yc  are  the  input  commands  and  0  and  7  are  the 

measured  pitch  and  yaw  angles,  respectively.  From  Equations 
(4.36)  and  (4.37),  the  differential  equations  for 
the  tracking  errors  are: 

£’e=-g  ;  (4.38) 

and 
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(4.39) 


Combining  Equations  (4.38)  and  (4.39)  with  the  time 
derivatives  of  Equations  ^4.34}  ana  (4.35),  a  state  equation 
that  is  appropriate  for  tracking  system  design  is  *£=4  H+jg  K, 

where  the  state  is: 


kt-  [£9  £9  q  f  6.  5e  5.  jr  (4.40) 

the  control  matrix  is  i;r=[u9  ur]7,  and  the  system  matrices 
are: 
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0 

-7  4  . 


(4.41) 


and 


r=00000000  2745 . 8 
“  ~tQQ0QQ00C  0 

In  a  manner  similar  to  that  used  in  the  roll  rate 
and  altitude  rate  controller  designs,  the  reduced  order 
model  is  formed  by  ignoring  the  servo  dynamics.  The 
resulting  state  is: 


0 

2745 


IT 


(4.42) 
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the  control  matrix  is  v;r=[6e  6,3  r,  and  the  system  matrices 
are: 

0  1  0  0  0  0 

0  0  -1  0  0  0 

00  0  00-1.62  (4.44) 

0  0  0  1  0  0 

0000-1  0 
0  0  1.78  0  0  0 

and 

_  [O  0  -11.19  0  0  0  T  (4.45) 

-  "  0  0  0  0  0  -12.28J 


2.  Pitch  and  Yaw  Angle  Controller  Design 

a.  Discretizing  the  Pitch  and  Yaw  Angle  System 
Using  the  sampling  period,  T,  equal  to  0.01 
seconds,  the  discrete-time  state  is: 

Kik)  T  =  [i?e(/c)  E9{k)  g(k)  E,[k)  £*  (k)  t(k)]T;  (4.46) 

the  control  matrix  is  v{k)  T  =  [5e(ic)  5,(ic)3T,  and  the 


discrete-time  system  matrices  are: 


Q(k)  =! 


1  0.01 

0  1 

0  0 

0  0 

0  0 

0  0 


-0 . 01 
0.9999 
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0.0178 
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0 

0101 

0 

0 
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n 

0 

0 

0.99 

0 


0  1 
0.0001  ! 
-0 . 0162 
0 
0 

0.9999 


(4.47) 


and 


ru)  = 


o  o 

0.0006  0 

■0.1119  0.001 

0  0 

0  0 

[-0.001  -0.1228J 


(4.48) 


b .  Gain  Determination 


The  LQR  weighting  matrices,  2  and  S'  for  the 
pitch  and  yaw  angle  control  loop  are  chosen  to  satisfy  the 
following  design  criteria: 

1.  The  overshoot  to  a  step  response  should  be  less  than  50 
percent . 

2.  The  five  percent  settling  time,  tj*,  is  less  than  or 
equal  to  2  seconds. 


The  allowable  overshoot  requirement  may  seem  too  liberal. 
However,  due  to  the  small  pitch  or  yaw  angles  (one  to  five 
degrees)  required  to  induce  translational  velocities  along 
the  body  fixed  y  and  z  axes,  the  overshoot  should  not 
present  any  problems  to  the  operator.  Computer  simulations 
have  indicated  that  excessive  control  values,  and  vane  angle 
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deflections  would  be  required,  if  the  controller  was 
designed  to  limit  the  overshoot  to  values  less  than  15 
percent.  The  large  overshoot  values  are  due  to  the  reduced 
order  controller.  If  a  full  order  controller  is 
implemented,  the  large  overshoot  values  could  be  eliminated 
with  small  vane  angles  and  a  reasonable  amount  of  control. 
However,  that  would  require  additional  sensors  and/or  an 
observer  as  discussed  previously. 

The  resulting  state  weighting  matrices  were  found 

to  be: 


6  0  0  0  0  0 

0  1  0  0  0  0 

00. 05  00  0 

*  0  0  0  6  0  0 

0  0  0  0  1  0 

00  0  0  0  .05 

The  steady-state  feedback  gain  matrix,  K,  is: 

K=  2.2741  1.8933  -0.6059  0.6390  0.5320  0.0024  ]  (4.50) 

“  [“0 .6390  -0.5320  -0.0024  2.2741  1.8933  -0.6059]' 

Figure  4.10  shows  the  MIMO  pitch  angle  and  yaw  angle 
controller  block  diagram  [Ref.  14:  p.  179]. 
c.  Simulation  Results 

Figure  4.11(a)  show  the  response  of  the  system  to 
a  commanded  pitch  angle  of  ten  degrees  and  a  commanded  yaw 
angle  of  zero.  The  large  value  of  overshoot  was  expected 
due  to  the  severe  coupling  of  the  axes.  Figure  4.11(b) 


S 


1  0 
0  1 


(4.49) 
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shows  the  induced  yaw  angle  from  the  commanded  pitch  angle. 
Figure  4.12  gives  a  similar  representation.  The  yaw  angle 
is  commanded  to  five  degrees  while  the  pitch  angle  is  held 
at  zero.  Finally,  Figure  4.13  shows  the  response  when  the 
pitch  and  yaw  angle  are  commanded  simultaneously  to  five 
degrees.  The  maximum  vane  deflection  applied  in  the 
simulations  documented  by  Figures  4.11,  4.12,  and  4.13  was 
1.7  degrees.  The  design  goals  for  overshoot  and  settling 
time  are  satisfied  for  the  linear  model  with  reasonable  vane 
deflections. 

d.  Singular  Value  Analysis 

Evaluating  the  robustness  of  the  roll  and 
altitude  rate  controllers  was  accomplished  by  computing 
their  respective  gain  and  phase  margins.  In  the  case  of  the 
MIMO  pitch  and  yaw  control  system,  the  description  of 
robustness  in  terms  of  the  gain  and  phase  margins  becomes 
more  complex.  Therefore,  a  different  approach  is  taken  in 
order  to  determine  the  robustness  of  the  pitch  and  yaw  angle 
control  system.  The  method  of  choice  is  to  apply  singular 
value  analysis  to  the  system  with  perturbations. 

Maciejowski  [Ref.  19]  presents  a  detailed  development  of  the 
theory  and  procedures  for  computing  the  singular  values  of  a 
MIMO  system  and  using  them  for  robustness  analysis. 
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Singular  value  analysis  is  used  to  examine  the 
robustness  of  the  MIMO  pitch  and  yaw  angle  controller.  Burl 
[Ref.  13]  demonstrates  that  a  system  is  stable  for  all  input 
multiplicative  perturbations,  a,  such  that 


H.  * 


1 

max  a  £(ju>) 

(J 


(4.51) 


where  o  is  the  largest  singular  value  of  the  matrix,  and 
T  ( jcu)  is  the  transfer  function  from  the  output  of  the 
perturbation  to  the  input  to  the  perturbation  as  shown  in 
Figure  4.14.  The  singular  values  for  the  MIMO  pitch  and  yaw 
angle  system  are  plotted  in  Figure  4.15.  The  largest 
singular  value,  a ,  is  equal  to  1.5.  Therefore,  applying 
Equation  (4.51),  the  MIMO  system  is  expected  to  be  robust 
for  perturbations  such  that  the  infinity  norm  of  a  is  less 
than  or  equal  to  0.667.  This  indicates  the  controller 
design  should  provide  good  results  when  applied  to  the 
nonlinear  model. [Pef.  13] 


E.  RESULTS  WITH  THE  NONLINEAR  SYSTEM 

The  controllers  designed  in  the  previous  sections  were 
designed  using  optimal  control  theory.  Because  the  servo 
states  were  ignored  in  determining  the  gains,  the  steady- 
state  gains  are  sub-optimal  in  respect  to  the  full  order 
system  models.  Through  simulations,  and  performance  and 
robustness  evaluations  it  was  determined  that  ignoring  the 
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servo  dynamics  would  have  no  adverse  affects.  This  section 
will  apply  the  controllers  designed  with  respect  to  the 
linear  system  models  to  the  nonlinear  model.  The  MATLAB 
programs  used  to  simulate  the  nonlinear  system  are  included 
as  Appendix  A. 

1.  Simulation  One  -  Figure  4.16 

The  system  was  commanded  to  climb  at  a  rate  of  ten 
feet/second  for  five  seconds.  The  slope  of  the  altitude 
plot  indicates  that  the  rate  of  ten  feet/second  is  achieved. 
The  roll  rate,  pitch  angle  and  yaw  angle  were  set  equal  to 
zero.  The  control  system  is  tracking  the  commanded  altitude 
rate  input  while  regulating  the  other  inputs  to  zero. 

Figure  4.16  shows  that  the  vehicle  climbs  to  an  altitude  of 
50  feet.  Note  the  coupling  between  the  roll  axis  (both  roll 
rate  and  roll  angle)  as  the  engine  accelerates.  The 
displacement  of  the  aileron  control  vanes  do  not  exceed  plus 
or  minus  one  degree  in  deflection. 

2.  Simulation  Two  -  Figure  4.17 

Figure  4 . 17  shows  the  response  to  a  commanded 
altitude  rate  of  25  feet/second  for  four  seconds.  A 
commanded  roll  rate  of  ten  degrees /second  for  two  seconds  is 
applied  starting  at  time  equal  to  six  seconds.  The  slope  of 
the  altitude  plot  indicates  that  the  altitude  rate  of  25 
feet/second  is  achieved.  The  vehicle  overshoots  t*e  desired 
100  feet  altitude  and  begins  to  approach  a  steady-state 
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Figure  4.16  Nonlinear  Simulation  One 
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value  of  approximately  80  feet.  The  large  overshoot  for  the 
altitude  is  due  to  the  momentum  of  the  vehicle.  Climbing  at 
a  high  rate  of  climb,  the  vehicle  is  expected  to  exhibit 
this  type  of  behavior.  Additionally,  the  control  system  is 
controlling  the  altitude  rate,  not  the  altitude.  The 
operator  would  have  to  account  for  this  in  his  control 
adjustments  during  flight.  Again,  the  coupling  between  the 
roll  axis  and  the  engine  acceleration  is  observed.  The 
commanded  roll  rate  drives  the  vehicle  to  an  roll  angle  of 
18°.  The  two  degree  error  is  due  to  the  minus  two  degrees 
of  roll  angle  caused  by  the  reactive  torque  to  the  altitude 
rate.  The  increased  aileron  vane  displacement  in  order  to 
perform  the  desired  maneuver  is  shown. 

3.  Simulation  Thraa  -  Figure  4.18 

The  vehicle  was  commanded  to  climb  at  a  rate  of  50 
feet/second  for  four  seconds.  A  pitch  angle  of  minus  five 
degrees  was  commanded  at  time  equal  to  six  seconds.  Figure 
4.18  illustrates  that  the  increased  rate  of  climb  increases 
the  altitude  overshoot.  The  coupling  between  the  engine  and 
roll  axis  is  present.  Additionally,  the  coupling  between 
the  pitch  and  yaw  axis  is  shown.  As  the  vehicle  is  pitched, 
the  yaw  angle  is  perturbed.  Note,  the  controller  is 
limiting  the  coupling  between  the  pitch  and  yaw  axis.  The 
yaw  angle  is  perturbed  less  than  one  degree  for  a  commanded 
pitch  angle  of  minus  five  degrees. 
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4.  Simulation  Four  -  Figure  4.19 

Figure  4.19  shows  the  system  response  for  an 
altitude  rate  of  50  feet/second  for  four  seconds.  The 
commanded  roll  rate  is  ten  degrees /second  for  one  second  at 
time  equal  to  three  seconds.  The  commanded  pitch  angle  is 
minus  three  degrees  applied  at  time  equal  to  six  seconds. 

The  control  vanes  displacement  angles  all  remain  less  than 
four  degrees  for  this  series  of  maneuvers.  The  controller 
is  able  to  provide  satisfactory  results  for  the  nonlinear 
model. 

F.  CONCLUSION 

The  controller  has  demonstrated  the  ability  to  control 
the  nonlinear  system  with  varying  inputs  applied.  It  can  be 
concluded  that  the  robustness  of  each  control  system  is 
sufficient  in  order  to  account  for  the  difference  in  system 
parameters  between  the  linear  models  and  the  nonlinear  model 
of  Archytas.  The  fact  that  the  control  systems  maintain 
stability  while  directing  commanded  inputs  to  steer  and 
maneuver  the  vehicle  supports  the  linear  modeling  approach 
taken  in  this  thesis.  The  next  step  in  the  design  process 
is  to  test  the  control  systems  on  a  prototype  vehicle  in  a 
controlled  experiment.  Such  a  test  was  performed  with  the 
roll  rate  control  system  and  an  Archytas  prototype  mounted 
on  a  test  stand.  The  results  of  this  test  will  be  discussed 
in  Chapter  V. 
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V.  CONCLUSIONS 


A.  ROLL  RATE  CONTROL  SYSTEM  FIELD  TEST 

Figure  5.1  shows  the  Archytas  prototype  (minus  the  wings 
and  canard)  mounted  on  a  test  stand  that  allows  the  vehicle 
to  spin  about  the  longitudinal,  or  x  axis.  This 
configuration  was  used  during  the  testing  of  the  roll  rate 
controller.  The  testing  was  accomplished  in  a  qualitative 
manner;  no  empirical  data  was  co'.lected,  the  object  of  the 
test  being  to  validate  the  roll  rate  controller  design.  The 
test  consisted  of  three  parts: 

1.  The  ability  of  the  roll  rate  controller  to  eliminate 
the  rotational  velocity  imparted  to  Archytas  from  the 
reactive  torques  applied  to  the  roll  axis  as  the  engine 
speed  is  varied  (decoupling) . 

2.  The  ability  of  the  roll  rate  controller  to  eliminate 
the  rotational  velocity  imparted  to  Archytas  from  the 
effect  of  cross-winds  (disturbance  rejection) . 

3.  The  ability  of  the  roll  rate  controller  to  allow  the 
operator  to  impart  or  terminate  a  rotational  velocity 
about  the  x  axis  (tracking  commanded  inputs) . 


Part  one  was  accomplished  by  varying  the  engine  speed 
from  idle  to  maximum  rpm  with  the  commanded  roll  rate  equal 
to  zero.  Each  time  the  engine  speed  was  varied,  the 
rotational  velocity  of  the  vehicle  was  eliminated  by  the 
roll  rate  controller.  It  was  observed  that  the  rotational 
velocity  was  eliminated  within  approximately  two  seconds. 
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Figure  5.1  Archytas  Prototype  Mounted 
on  the  Test  Stand 
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This  response  is  consistent  with  both  the  linear  and 
nonlinear  simulations. 

Part  two  was  accomplished  by  setting  the  engine  speed  at 
constant  values  ranging  from  idle  to  maximum  rpm  with  the 
commanded  roll  rate  equal  to  zero.  A  disturbance  about  the 
x  axis  was  imparted  to  Archytas  by  pushing  the  test  stand 
mounting  bracket  connected  to  the  vehicle.  The  effect  of 
pushing  the  connecting  bracket  is  analogous  to  a  rotational 
velocity  imparted  by  cross-winds.  It  was  observed  that  the 
. jtational  velocity  was  eliminated  within  approximately  two 
seconds.  This  response  is  consistent  with  both  the  linear 
and  nonlinear  simulations. 

Part  three  was  accomplished  by  inputting  a  commanded 
roll  rate  while  the  engine  speed  was  held  constant  and  also 
varied.  It  was  observed  that  the  operator  could  position 
Archytas  at  a  desired  angle  under  either  test  condition. 
Additionally,  disturbances  were  imparted  to  Archytas  during 
these  tests.  The  roll  rate  controller  proved  robust  enough 
to  eliminate  the  disturbances  and  allow  the  operator  to 
position  the  vehicle.  This  response  is  consistent  with  both 
the  linear  and  nonlinear  simulations. 

Based  on  the  computer  simulations,  both  linear  and 
nonlinear,  and  the  results  of  the  above  tests,  it  is 
concluded  that  the  roll  rate  controller  design  is  valid. 
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B.  FUTURE  RESEARCH 

The  development  of  the  Archytas  control  systems  was 
accomplished  using  the  data  from  the  AROD  as  first 
approximations.  Future  research  will  include  the 
computation  of  Archytas  specific  data  through  empirical 
methods  or  obtaining  the  data  by  wind  tunnel  testing.  This 
would  allow  for  a  better  model  of  the  Archytas  to  be 
developed. 

The  second  order  dynamical  model  of  the  servos  was 
obtained  with  the  vanes  mounted  on  the  servos  with  zero 
downwash  from  the  propeller.  Future  research  should  include 
modeling  of  the  servos  with  the  engine  at  hover  rpm  and  the 
vanes  positioned  within  the  downwash. 

The  effects  of  inclusion  of  the  vane  position  angles  in 
the  control  systems  could  be  evaluated  further.  In 
particular,  the  effect  the  vane  angle  position  feedback 
would  have  in  the  reduction  of  the  overshoot  in  the  pitch 
and  yaw  angle  controller  should  be  investigated. 
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APPENDIX  A 


MATLAB  SIMULATION  PROGRAMS 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  PROGRAM  NAME:  nonlin.m 

% 

%  DESCRIPTION: 

%  This  program  will  simulate  the  Archytas  nonlinear  model% 
%  with  the  linear  controllers.  The  Control  Laws  are  the  % 
%  Steady-state  Linear  Quadratic  Regulator  solutions  with  % 
%  the  necessary  modifications  to  the  performance  measure. % 
%%%%%%%%%%%*,%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
clear 

%  Call  the  program  noninp.m  to  enter  the  desired  roll  rate, 

%  alttitude  rate,  pitch  and  yaw  angles. 

noninp 

%  Call  the  program  in_cond.m  to  load  the  initial  conditions 

and  %  constants  for  the  simulation. 

incond 

%  Call  the  program  n_l_inp.m  to  load  the  wind  tunnel  data. 
n_l_inp 

%  Loop  for  simulation 
for  k*l :num_of  steps 1 

%  Compute  the  speed  of  the  propeller 
speed (k)  =  speed_hover  +  del_rpm(k) ; 

%  Compute  the  thrust  supplied  by  the  engine  to  the  vehicle 
thrust  (k)  =  thrust_hover  +  xrpm  *  del_rpm(k) ; 

%  Check  thrust  limits 
if  thrust (k)  <  thrust_min 
thrust (k)  =  thrust_min; 
elseif  thrust (k)  >  thrust_max 
thrust (k)  =  thrustjnax; 

end 

velocity_tip(k)  =  speed  hover  +  del_rpm(k) ; 
del_tip(k)  =  velocity_tIp(k)  -  speed_hover; 

tip_del(k)  =  velocity _tip(k)  -  speed(k); 

stad(k)  =  tip_del(k); 

%  Compute  vehicle  velocity  and  angle  of  attack  total 
aoatot(k)  *  alphajraax; 

velocity^tot (k) =sqrt(u(k) ^2+v(k) A2+w(k)  A2) ; 
if  velocTty_tot (k)  >  vaero 

vwterm  =  sqrt(v(k)A2  +  w(k)A2); 
aoatot(k)  *  asin(vwterm/velocity_tot (k) ) ; 
if  aoatot(k)  <  alpha_min 
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aoatot(k)  *  alphajnin; 
elseif  aoatot(k)  >  alphajmax 
aoatot(k)  =  alpha_max; 

end 

end 

veq(k)  *  tablel(veq  table, aoatot(k) ) ; 
newveq(k)  =  sqrt (weIght_ratio) *veq(k) ; 
rhova(k)  *  rho  *  newveq(k)  *  pi; 
velocity_tip_2 (k)  =  2.0/velocity_tip(k) ; 
velocity_delta(k)  -  velocity_tot (k)  -  newveq(k) ; 

%  This  section  assumes  that  no  aerodynamic  forces  and  moments 
%  exist 

if  velocity _tot(k)  <=  vaero 

fax  =  thrust (k)  *  gravity /weight ; 

fay  =  0.0; 

faz  *  0.0; 

mrx  =  0.0; 

mpy  =  0.0; 

myz  =  0.0; 

pitch_trim  »  0.0; 

yaw_trim  =  0.0; 

pitch_factr  =  1.0; 

yaw_factr  =  1.0; 

elseif  velocity _tot (k)  >  vaero 
if  aoatot(k)  <=  alpha_min 

fax  =*  thrust (k)  *  gravity /weight; 

fay  »  0.0; 

faz  =  0.0; 

mrx  =  0.0; 

mpy  *  0.0; 

myz  =  0.0; 

pitch_trim  =  0.0; 

yaw_trim  *  0.0; 

pitch_factr  =  l.o; 

yaw_factr  =  1.0; 

end 

else 

%  This  section  computes  the  aerodynamic  forces  and  moments 

lleq  =  weight_ratio  *  leq; 

ddeq  *  weight~ratio  *  deq; 

sseq  =  weight_ratio  *  seqq; 

req  =  tablel(req_table,aoatot(k) ) ; 
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peq  =*  tablel (peq_table, aoatot (k) ) ; 
yeq  *  tablel (yeq_table, aoatot (k) ) 

rreq  »  weight_ratio  *  req; 

ppeq  =  weigth_ratio  *  peq; 

yyeq  =>  weight_ratio  *  yeq; 

cldel  =  tablel(cldeX_table,aoatot(k) ) ; 
cddel  =*  tablel (cddel__table, aoatot (k) )  ; 
csdel  *  tablel (csdel_table, aoatot (k) ) ; 

Xdel =  velocity_tip_2 (k) *lleq-rhova*weight_ratio*cldel; 
ddel  =  velocity__tip_2  (k)  *ddeq-rhova*weight_ratio*cddel; 
sdel = velocity~tip_2 (k) *sseq-rhova*weight_ratio*csdel ; 

crdel  *  tablel (crdel_table, aoatot (k) ) ; 
cpdel  «  tablel (cpdel_table, aoatot (k) ) ; 
cydel  =  tablel (cydel_table, aoatot (k) ) ; 

rdel  =  velocity__tip_2  <k) *rreq-rhova*weight_ratio*crdel ; 
pdel  »  velocity ”tip_2 (k) *ppeq-rhova*weight_ratio*cpdel; 
ydel  =  velocity_tip_2 (k) *yyeq-rhova*weight_ratio*cydel; 

Xslope  =  tablel (lslope_tablef aoatot (k) ) ; 
dslope  =  tablel (dslope_table, aoatot (k) ) ; 
sslope  *  tablel (sslope_table, aoatot (k) ) ; 

ffl  -  lleq+ldel*del_tip(k) +lslope*weight_ratio* 
velocity_delta(k) ; 

ffd  =  ddeq+ddel*del_tip(k)+dslope*weight_ratio* 
velocity_delta(k) ; 

ffs  *  sseq+sdel*del_tip(k) +sslope*weight_ratio* 
velocity_delta(k) ; 

rslope  =  tablel (rslope_table, aoatot (k) ) ; 
pslope  *  tablel (pslope_table, aoatot (k) ) ; 
yslope  =  tablel (yslope_table, aoatot (k) ) ; 

mr  =  rreq+rdel*del_tip(k)+rslope*weight_ratio* 
velocity_delta(k) ; 

mp  =  ppeq+pdel*del_tip(k)+pslope*weight_ratio* 
velocity_delta(k) ; 

my  =  yyeq+ydel*del_tip(k) +yslope*weight_ratio* 
velocity_delta(k) ; 

fs  *  ffs  *  gravity/ weight ; 

fl  =  ffl  *  gravity /weight; 

fd  *  ffd  *  gravity /weight; 

delta  *  atan(v(k) /w(k) ) ; 

fax  *  fl  *  sin(aoatot (k) )  -  fd  *  cos(aoatot (k) ) ; 


101 


fayz  =-fl  *  cos(aoatot(k) )  -  fd  *  sin(aoatot (k) ) ; 
fay  *  fayz  *  sin(delta) ; 

faz  =  fayz  *  cos (delta); 

mrx  -  my  *  sin(aoatot (k) )  -  mr  *  cos (aoatot (k) ) ; 

mpyyz  =-my  *  cos (aoatot (k) )  -  mr  *  sin(aoatot (k) ) ; 
mpy  =  mpyyz  *  sin(delta) ; 

myz  =  mpyyz  *  cos(delta); 

velocity_trim  = tablel(velocity_trim_table, aoatot (kl) ) ; 
vaneff  *  tatlel ( vaneff _table, aoatot (k) ) ; 


pitch_trim  =  velocity_trim  *  cos(delta); 
yaw_trim  =-velocity_trim  *  sin (delta) ; 
pitch_factr  =  vaneff;  ~ 
yaw_factr  =  vaneff; 


end 

lat(k)  =  ( (rscale*mrx*gravity*144 . 0) /Ixx) +(La*del_a(k) ) - 
(velocity_delta(k) *2  /Ixx)  *  proptorq; 

mat(k)  «  ( (mpy*gravity*144 . 0) /Iyy) +(Me*pitch_factr* 

(del_e (k) -pitch_trim) ) ; 

nat(k)  =  ( (myz*gravity*144.0) /Izz) +(Nr*yaw_factr  * 

(del_r (k) -yaw_trim) ) ; 

%  Begin  Differential  Equations  of  Motion 

%  Altitude  Differential  Equations 
tempi (k)*u(k) *cos (theta (k) ) *cos(psi(k) ) ; 
temp2 (k)=v(k) * (sin(phi (k) ) *sin(theta(k) ) *cos (psi (k) ) - 
cos (phi (k) ) *sin(psi (k) ) ) ; 

temp3 (k)=w(k) * (cos (phi (k) ) *sin (theta (k) ) *cos(psi(k) ) + 
sin(phi (k) ) *sin(psi (k) ) ) ; 
u_earth (k) =  tempi (k) +temp2 (k) +temp3 (k) ; 

%  Heading  Angle  Differential  Equations 
temps (k)=p(k) *cos (theta (k) ) *cos (psi (k) ) ; 
temp6(k)*q(k) * (sin (phi (k) ) *sin (theta (k) ) *cos (psi (k) ) - 
cos (phi (k) ) *sin(psi(k) ) ) ; 

temp7(k)*r(k)*(cos(phi(k) ) *sin(theta (k) ) *cos(psi(k) ) 
+sin(phi(k) ) *sin(psi(k) ) ) ; 
hdg_dot (k) ®temp5 (k) +temp6 (k) +temp7 (k) ; 

%  Euler  Angle  Differential  Equations:  Pitch,  Yaw,  and  Roll 
theta_dot(k)=q(k) *cos(phi(k) )-r(k) *sin(phi(k) ) ; 
phi_dot=p(k) +(q(k) *sin(phi(k) )+r(k)*cos(psi(k) ) )* 
tan (theta (k) ) ; 
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psi_dot(k)-(q(k)*sin(phi(k) )+r(k) *cos(phi(k) ) ) / 
cos (theta (k) ) ; 

%  Vehicle  Pitch,  Yaw  and  Roll  Rate  Differential  Equations  WRT 
%  Body  Axes 

p_dot(k)  =  (bca*q(k) *r (k) ) -(Iralph*stad(k) )+lat(k) ; 
q_dot(k)  =  (cab*p(k) *r (k) ) -(Irx*r (k) *speed(k) /Iyy) +mat (k) ; 
r_dot(k)  -  (abc*p(k) *q(k) ) +(Irx*q(k) *speed(k) /Izz) +nat (k) ; 

%  Vehicle  Velocity  Differential  Equations  WRT  Body  Axes 
u_dot(k)=( (v(k) *r(k) ) -(w(k) *q(k) ) ) -gravity*cos (theta (k) ) * 
cos(psi(k) )+fax; 

v_dot (k)=( (w(k) *p(k) ) -(u(k) *r (k) ) ) +gravity*cos (theta (k) ) * 
sin(psi(k) )+fay; 

w_dot(k)=( (u(k) *q(k) )  -(v(k) *p(k) ) ) -gravity* 
sin (theta (k) )+faz; 


%  Define  Error  Variable  in  Pitch  and  Yaw 
E_theta(k)=theta_com(k) -(theta (k) +noise(k) ) ; 

E_psi (k) =psi_com(k) -(psi (k) +noise(k) ) ; 

%  Define  Error  Variable  in  Roll 
E_P(k)=P_coin(k)‘"(P(J<:)+noise(^) ) ! 

%  Define  Error  Variable  in  Altitude 
E_h_dot  (k)  =h_dot_com(k)  -  (u(k)  +noise(k) )  ; 

%  Integrate  Errors 

E_thetain (k+l) =E_thetain (k) +Ts*E_theta (k) ; 

E_psiin(k+l)=E  psiin(k)+Ts*E_psi(k) ; 
E_pin(k+l)=E_pin(k)+Ts*E_p(k) ; 

E_h_dotin(k+l) =E_h_dotin(k) +T s *  E_h_dot ( k ) ; 

%  Define  Control  System  Gains 
%  Roll  Rate  Control  Gains 
Kl=[  0.4376  -0.2027]; 

%  Pitch  and  Yaw  Angle  Control  Gains 

K2=(  2.2741  1.8933  -0.6059  0.6390  0.5320  0.0024; 

-0.6390  -0.5320  -0.0024  2.2741  1.8933  -0.6059]; 

%  Altitude  Rate  Control  Gains 
K3=[ -0 . 3060  0.2003  0.0038]; 

%  Calculate  the  Aileron  (roll)  Servo  Control  Input 
Ua(k)='-Kl*[E_pin(k+l) ;  (p(k) +noise(k) )  ] ; 

%  Calculate  the  Elevator  (pitch)  Servo  Control  Input 
Ue(k)=-K2 (1, : ) *[E_thetain(k+l) ;E_theta(k) ;q(k) ;E_psiin(k+l) ; 

E_psi(k) ;r(k) ] ; 
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%  Calculate  the  Rudder  (yaw)  Servo  Control  Input 
Ur(k)=-K2(2, : )  *[E_thetain(k+l)  ;E_theta(k)  ;q(k)  ;E_psiin(k+l) 

E_psi(k) ;r(k) ] ; 

%  Calculate  the  Throttle  (altitude)  Servo  Control  Input 
Ut (k)=-K3*[E_h_dotin(k+l) ;u(k) ;del_rpm(k) ]; 

%  Begin  Integration  of  Equations  of  Motion 
del_a_dot_dot (k) =-Hl*del_a_dot (k) -H2*del_a (k) +H2*Ua (k) ; 
del  e_dot_dot(k)=-Hl*del  e_dot(k)-H2*del_e(k)+H2*Ue(k) ; 
del_r_dot_dot (k) =-Hl*del~r_dot (k) -H2*del_r (k) +H2*Ur (k) ; 
de l_t_do t_dot (k) =-Hl*del_t_dot (k) -H2*del_t (k) +H2*Ut (k) ; 

%Roll,  pitch,  yaw 
p(k+l)=p(k)+Ts*p_dot(k) ; 
q(k+l)=q(k) +Ts*q_dot(k) ; 
r (k+l)=r (k)+Ts*r_dot(k) ; 

%Velocities 

u(k+l)=su(k)  +Ts*u_dot(k) ; 
v(k+l) =v(k) +Ts*v_dot (k) ; 
w(k+l)=w(k)+Ts*w_dot(k) ; 

%Euler  angles 

theta (k+1) =theta (k) +Ts*theta_dot (k) ; 
phi(k+l)=phi(k)+Ts*phi_dot(k) ; 
psi (k+1) =psi (k) +Ts*psi_dot (k) ; 

%Servos 

del_a_dot (k+1) =del_a_dot (k) +Ts*del_a_dot_dot (k) ; 
if  del_a_dot (k+1)  >  rmax 
de l_a_dot ( k+ 1 )  =  rmax ; 
elseif  del_a_dot (k+1)  <  -rmax 
del_a_dot(k+l)  =  -rmax; 

end 

del_e_dot(k+l)=del_e_dot(k) +Ts*del_e_dot_dot (k) ; 
if  del_e_dot (k+1)  >  rmax 
del_e_dot (k+1)  «  rmax; 
elseif" del_e_dot (k+1)  <  -rmax 
del_e_dot (k+1)  =  -rmax; 

end 

del_r_dot (k+1) =del_r_dot (k) +Ts*del_r_dot_dot (k) ; 
if  del_r_dot (k+1)  >  rmax 
del_r_dot (k+1)  =  rmax; 
elseif  del_r_dot (k+1)  <  -rmax 
de l_r_do t ( k+ 1 )  =  -rmax; 

end 
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del_t_dot (k+l) *del_t_dot (k) +Ts*del_t_dot_dot (k) ; 
if  del_t_dot (k+l)  >  rmaxt 
del_t_dot (k+l)  =  rmaxt; 
elseif  del_t_dot(k+l)  <  -rmaxt 
del_t_dot (k+l)  =  -rmaxt; 

end 

%  Engine 

del_rpm_dot (k+l) =-w_e*del_rpm(k) +Ke*w_e*del_t (k) ; 

del_a (k+l) =del_a (k) +Ts*del_a_dot (k); 

del_e(k+l)=del_e(k)+Ts*del_e_dot(k) ; 
if  del_e(k+l)  >  maxdfl 
de l_e ( k+ 1 )  *  maxdf 1 ; 
elseif  del_e (k+l)  <  -maxdfl 
del_e(k+l)  =  -maxdfl; 

end 

del_r (k+l)-del_r (k)+Ts*del_r  dot(k) ; 
if  del__r  (k+l)  >  maxdfl 
del_r(k+l)  *  maxdfl; 
elseif  del_r(k+l)  <  -maxdfl 
del_r(k+l)  =  -maxdfl; 

end 

del_t  (k+l)  =del_t  (k)  +Ts*del_t_dot  (k) ; 
if  del_t(k+l)  >  tamax 
del__t(k+l)  *  taunax; 
elseif  del_t(k+l)  <  -tamax 
del_t(k+l)  *  -tamax; 

end 


%  Engine 

del_rpm(k+l)  =*del_rpm(k)  +Ts*del_rpm_dot (k)  ; 
%  Altitude 

alt (k+l) =alt (k) +Ts*u_earth (k) ; 

%  Heading  Angle 

heading (k+l) =heading(k) +Ts*hdg_dot (k) ; 

t(k)=Ts*(k-l) ; 

%  End  of  simulation  loop 
end 

t (k+l) =Ts*k; 

u_earth(k+l)*u_earth(k) ; 
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%  Plot  the  Results  of  the  Simulation 

subplot  (  221)  ,  plot  (t  ,  alt)  ,grid,xlabel  (  'Time 
(sec) ' ) ,ylabel( 'Altitude  (ft) ') ; 
title ( 'Altitude' ) ; 

subplot ( 222 ) , plot ( t , p*57 .2958) , grid, xlabel ( 'Time 
(sec) ') ,ylabel( 'Roll  Rate  (deg/sec)'); 
title ('Roll  Rate  -  p'); 

subp 1 o t (  2  2  3  )  ,plot(t,phi*57 . 2958)  ,grid,xlabel('Time 
(sec) ') , y label (' Angle  (deg)'); 
title( 'Roll  Angle  Phi'); 

subplot (224) ,plot(t,theta*57 .2958) , grid, xlabel ( 'Time 

(sec) ') ,ylabel( 'Angle  (deg)'); 

title ('Pitch  Angle  Theta'); 

meta  nonl 

subplot (ill) ; 

subp lot (221) , plot (t , psi*57 . 2958), grid, xlabel ( 'Time 
(sec) ') ,ylabel ( 'Angle  (deg)'); 
title ('Yaw  Angle  Psi'); 

subplot (222) ,plot(t,del_a*5  7.2  9  58)  ,  grid, xlabel ( 'Time 
(sec) ') ,ylabel( 'Angle  (deg)'); 
title ( 'Aileron  Vane  Angle'); 

subplot (223) , plot ( t , del_r*57 .2958)  ,grid,xlabel( 'Time 
(sec) ') ,ylabel ( 'Angle  (deg)'); 
title ( 'Rudder  Vane  Angle'); 

subplot (224) , plot (t,del_e* 57 . 2958) , gr id , xlabe 1 ( 'Time 

(sec) ') ,ylabel( 'Angle  (deg)'); 

title ( 'Elevator  Vane  Angle'); 

meta  non2 

subplot (111) ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  PROGRAM  NAME:  noninp.m  % 
%  % 
%  DESCRIPTION:  % 
%  This  program  prompts  the  user  for  the  desired  % 
%  simulation  inputs.  % 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Prompt  the  user  for  the  desired  simulation  time. 
disp{'  Enter  the  simulation  time  in  seconds.'); 
simulation_time  =  input ('>»>  '); 

%  Prompt  the  user  to  enter  to  desired  sampling  time  -  Ts. 
disp('  Enter  the  desired  sampling  time  -  Ts.'); 
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Ts  -  input  {'»»  '); 

%  Compute  the  number  of  simulation  steps  based  on  the 
desired  simulation  time. 

num_of_stepsl  =  round (simulation_time/Ts) ; 

%  Prompt  the  user  for  the  length  of  the  altitude  rate  step 
input  time. 

disp('  Enter  the  length  of  the  altitude  rate  input  in 
seconds . ' ) ; 

step_input_timel  =  input  ('»>>  '); 

%  If  step_input  entered  is  greater  than  simulation_time 
entered,  set  them  equal, 
if  step_input_timel  >  simulation_time 
step_input_timel  =  simulation_time; 

end 

if  step_input_timel  ~=  0 

%  Prompt  the  user  for  the  desired  magnitude  of  the  step 
%  input  in  #  of  deg/sec. 

disp('  Enter  the  magnitude  of  the  altitude  rate  input  in 
feet /sec. ' ) 

h  dot  =  input  ( ' »»  ' ) ; 
elseTf  step_input_timel  ==  0 
h_dot  =0.0; 

end 

%  Compute  the  number  of  step  input  steps  based  on  the 

%  desired  step  input  time. 

num_o  f _s t eps  2  =  step_input_timel/Ts; 

num_of_steps2  =  round (num_of_steps2) ; 

%  Compute  the  number  of  zeros  to  be  padded  to  the  input 
%  based  on  the  length  simulation  time  and  length  of  the 
%  desired  step  input  time. 

num__of_zerosl  =  num_of_stepsl  -  num_of _steps2 ; 

%  Define  the  commanded  input  vector 
h_dot_com  =  [h_dot  *  ones  (1,  num__of__steps2) 
zeros ( 1, num_of_zerosl)  ] ; 

%  Prompt  the  user  for  whether  he  wants  to  include  a  roll 
rate  in  the  simulation. 

ans  =  input (' Do  you  want  to  include  a  roll  rate  in  the 
simulation?  y/n  [y] :  ','s'); 
if  isempty(ans) 
ans  ='y' ; 

end 

if  ans  =='y ' 
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%  Prompt  the  user  for  the  start  time  for  the  roll  rate 
step  input. 

disp('  Enter  the  start  time  for  the  roll  rate  input.'); 
start_timel  =  input ('»»  '); 

%  Prompt  the  user  for  the  length  of  the  roll  rate  input. 
disp('  Enter  the  length  of  the  roll  rate  imput  in 
seconds. ' ) ; 

length_roll_rate  =  input ('>>>>  '); 

%  Compute  the  number  of  zeros  to  pad  front  of  the 
%  command  vector. 
num_of_zeros2  =  start_timel/Ts ; 
num_of_zeros2  =  round (num_of_zeros2 ) ; 

%  Computer  the  number  of  step  input  steps  and  number  of 
%  zeros  to  pad  back  of  the  command  vector 

num_of_steps3  -  length_roll_rate/Ts; 
num_of_steps3  =  round(num_of_steps3) ; 
num_of_zeros3  =  num_of_stepsl- 

round( (start_timel+length_roll_rate) /Ts) ; 

%  Prompt  the  user  for  the  desired  magnitude  of  the  roll 
%  rate  input  in  #  deg/se. 

disp(' Enter  the  magnitude  of  the  roll  rate  in 
degrees/ sec. ' ) ; 
pi  =  input ('»»  '); 

%  Define  the  commanded  input  vector 

p_com  =  [ zeros (l,num_of_zeros2)  pi  *  0.017453  * 

ones ( 1 , num_of _steps3 )  zeros ( 1 , num_of_zeros3 ) ] ; 

else 

p_com  =  [ zeros ( l, num_of_stepsl) ] ; 

end 

%  Prompt  the  user  for  whether  he  wants  to  include  a  pitch 
%  angle  in  the  simulation. 

ans  =  input ('Do  you  want  to  include  a  pitch  angle  in  the 
simulation?  y/n  [yj :  ','s'); 
if  isempty(ans) 
ans  =  ' y ' ; 

end 

if  ans  -='y ' 

%  Prompt  the  user  for  the  start  time  for  the  pitch  angle 
%  command. 

disp('  Enter  the  time  at  which  you  desire  to  input  the 
pitch  angle.') 

start_time2  =  input ('»»  '); 

%  Compute  the  number  of  zeros  to  pad  front  of  the 
%  command  vector 

num_of_zeros4  =  round ( start_time2 /Ts ) ; 
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%  Compute  the  number  of  ones  to  multiply  the  desired 
%  angle  by  to  create  the  command  vector. 
number_of_onesl  =  num_of_stepsl  -  round (start_time2/Ts) ; 
%  Prompt  the  user  for  the  desired  pitch  angle  in 
%  degrees. 

disp('  Enter  the  desired  pitch  angle  in  degrees.'); 
theta_com  =  input ('>>>>  '); 

theta_com  =  [ zeros (l,num_of_zeros4)  theta_com  *  0.017453 
*  ones (1, number _of  onesl) ] ; 

else 

theta_com  =  [ zeros ( 1, num_of_stepsl) ] ; 

end 

%  Prompt  the  user  for  whether  he  wants  to  include  a  yaw 
%  angle  in  the  simulation. 

ans  =  input ('Do  you  want  to  include  a  yaw  angle  in  the 
simulation?  y/n  [y] :  ','s'); 
if  isempty(ans) 
ans  -  'y' ; 

end 

if  ans  == ' y ' 

%  Prompt  the  user  for  the  start  time  for  the  yaw  angle 
%  command . 

disp('  Enter  the  time  at  which  you  desire  to  input  the 
yaw  angle . ' ) 

start_time2  =  input ('»»  '); 

%  Compute  the  number  of  zeros  to  pad  front  of  the 
%  command  vector 

num_of_zeros4  -  round(start_time2/Ts) ; 

%  Compute  the  number  of  ones  to  multiply  the  desired 
%  angle  by  to  create  the  command  vector. 
number_of_onesl  =  num_of_stepsl  -  round (start_time2/Ts) ; 
%  Prompt  the  user  for  the  desired  pitch  angle  in 
%  degrees . 

disp('  Enter  the  desired  yaw  angle  in  degrees.'); 
psi_com  =  input ('»»  '); 

psi_com  -  ( zeros (l,num_of_zeros4)  psi_com  *  0.017453  * 
ones ( 1 , number_of _onesl) ] ; 

else 

ps i_com  =  [ zeros ( 1 , num_of _steps 1 ) ] ; 

end 

%  Prompt  the  user  for  whether  he  wants  to  add  sensor  noise 
%  to  the  simulation. 

ans  *  input  (/r>o  you  want  to  include  sensor  noise  in  the 
simulation?  y/n  [y] :  ','s'); 
if  isempty(ans) 
ans  =  'y' ; 

and 

if  ans  -«'y' 
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the  random  vector  to  represent  the  measurement 


%  Define 
%  noise. 

%  Assume  the  accuracy  is  +/-  i.745e-04 
rand ( 'uniform') ; 
noise  =  rand (1 , num_of_stepsl) ; 
vard  =  (1/12)  *  ( 1 . 745e-04 ) A  2 ;  % 

A  =  sqrt (vard/cov (noise) ) ;  % 


else 


noise  =*  A  .*  noise; 

%  Adjust  the  mean  to  zero 
noise  =  noise  -  mean (noise) 


end 


noise  =  [ zeros ( 1, num_of_stepsl) ) ; 


Desired  variance 
Scalar  multiplier  to 
change  the  variance 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  PROGRAM  NAME:  in  cond.m  % 
%  % 
%  DESCRIPTION:  % 
%  This  program  inputs  the  initial  conditions  and  % 
%  constant  values  used  in  nonlin.m  % 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  CONSTANTS 

vaero=0.5;  speed_max=859.0;  speed_min=!649 . 0 ; 
prop_torq=0.0729;  thrust_min=35.0;  thrust_max=105 . 0; 
tamax=100 . 0 ;  rmaxt=100. 0; 
maxdf 1=0. 5236;  rmax=0. 87266; 

La=-21.04;  Me=-9.01;  Nr=-ll.40; 

Hl-71.1;  H2=2745.8; 

W_e=2.0;  Ke=837 . 758 ;  Xrpm=0.2122;  speed_hover=712 . 0943 ; 
alpha_min=0 . 174533 ;  alpha_max«l . 570796 ; 
gravity*32.l74;  pi=3 . 1415962 ;  rad_to_deg=180. 0/pi; 
Ixx»6908.4;  Iyy-22944 . 64 ;  Izz=21515.64;  Irx=4 1.62; 
bca*0. 20685;  cab*o. 63663;  abc*-o. 74533;  lralph»0. 00916; 
weight=85.0;  thrust_hover=85 . 0;  velocity_tip=*722 . 5663 ; 
rho«0. 00192 ;  weight_ratio*l. 0; 

leq=85.0;  deq=0.0;  seqq=0.0;  sslope=0 . 0801 ;  _*scale=0.0; 

%  INITIAL  CONDITIONS 


p(l)=0.0; 
q ( 1 ) =0 . 0 ; 
r(l)*0.0; 
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u(l)=0.0; 
v ( 1) *0 . 0 ; 
w(i)*0.0; 

phi(l)  =0.0; 
psi(l)  =0.0; 
theta (1) =0.0; 

alt ( 1)  =0.0; 

dist(l)  =0.0; 
heading (1)  =0.0; 

E_thetain ( 1 )  =0.0; 

E_psiin(l)  =0.0; 

E_pin(l)  =0.0; 

E_h_dot in ( l )  =0.0; 

del_a(l)  =0.0; 
del_e(l)  =0.0; 
del_r(l)  =0.0; 
del_t (1)  =0.0; 

del_rpm(l)=0.0; 

del_a_dot ( l ) =o . 0 ; 
del_e_dot ( 1 ) =0 . 0 ; 
del_r_dot ( l ) =0 . 0 ; 
de l_t_dot ( l ) =0 . o ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  PROGRAM  NAME:  n_l_inp.ro  % 
%  % 
%  DESCRIPTION:  % 
%  This  file  inputs  the  wind  tunnel  data  into  the  % 
%  MATLAB  wordspace  % 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%«%%%%%%%%%%%«%%%%%%%%%%%%%%%% 

aoatotl=[0. 174533  0.349066  0.523599  0.785398  0.872665... 

0.959931  1.047198  1.134464  1.221730  1.308997... 
1.570796]; 

CRdell  =[0.00  0.00  0.03  0.06  0.06 

0.07  0.08  0.07  0.07  0.07  0.05]; 

Rslopel=[0. 00  0.01  0.07  0.14  0.16 

0.18  0.18  0.18  0.17  0,18  0.20]; 
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Reql  =[2.16  4.84  6.66  7.23  7.76 

3.66  8.92  6.89  6.66  5.68  0.00]; 

CPdell  =[-0.08  -0.09  -0.08  -0.02  -0.00 

0.03  0.07  0.12  0.08  0.12  0.00]; 

Pslopel=[-0.29  -0.24  -0.21  -0.04  -0.01 

0.08  0.17  0.29  0.19  0.31  0.40]; 

Peql  =[-5.06  2.24  7.46  15.32  16.65  17.12 
17.37  17.21  18.03  16.10  0.00]; 

CYdell  =[0.03  0.04  0.04  0.03  0.03 

0.03  0.03  0.01  0.01  0.01  0.00]; 

Yslopel=[0. 12  0.12  0.11  0.08  0.07 

0.06  0.06  0.02  0.02  0.02  0.00]; 

Yeql  =[4.54  6.36  6.82  3.13  3.05 

3.73  3.17  2.34  2.34  1.81  1.00]; 

VanEf f 1=[ 1 . 50  1.45  1.40  1.35  1.30 

1.25  1.20  1.15  1.10  1.05  1.00]; 

veql  =[172.4  115.4  87.65  62.85  55.45 

48.00  41.2234.68  28.98  22.37  0.0]; 


CLdell  =[0.16  0.28  0.32  0.29  0.30  0.26 

0.23  0.19  0.20  0.16  0.00]; 

Lslopel=[ 0 . 60  0.81  0.82  0.73  0.74  0.64 

0.56  0.46  0.50  0.40  0.00); 

COde 11  =[0.44  0.46  0.46  0.48  0.50  0.50 

0.50  0.50  0.54  0.52  0.50); 

Dslopel=[1.67  1.32  1.19  1.20  1.23  1.22 

1.24  1.23  1.34  1.32  1.25]; 

CSdell  =[-0.01  -0.02  -0.03  -0.01  -0.01 

-0.02  0.01  0.00  0.02  0.04  0.06]; 

Vtriml  =[0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0] 


crdeltable  =[aoatotl'CRdell' ] ; 
rslope_table=[ aoatot l'Rs lope 1' ] ; 
req_table  =[aoatotl'  Reql']# 


cpdel_table  =[aoatotl'  CPdell']? 
pslope_table=[aoatotl'  Pslopel' ] ; 
peq_table  =[ aoatot 1'  Peql']; 

cydel_table  =[aoatotl'  CYdell']; 
yslope_table=[aoatotl'  Yslopel' ] ; 
yeq_table  = [ aoatot 1 '  Yeql ' ] ; 


vanEf f_table=[aoatotl'  VanEf f 1' ] ; 
veq_table  =[ aoatot 1'  Veql']; 
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cldel_table  =[aoatotl'  CLdell'}; 
lslope_table= [ aoatotl '  Lslopel ' ] ; 

cddel_table  =[aoatotl'  CDdell'J; 
ds lope_table= [ aoatot l '  Dslopel' ] ; 

csdel^table  =[ aoatotl'  CSdell']; 
velocTty_trim_table  -[aoatotl'  Vtriml']; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
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APPENDIX  B  CONTROL  SERVOS 
The  control  servos  incorporated  into  Archytas  are 
position-commanded,  quarter-scale  model  airplane  servos. 

Five  of  these  servos  are  used  on  Archytas:  four  to  drive  the 
aerodynamic  control  surfaces  and  one  to  drive  the  engine 
throttle.  These  servos  can  be  modeled  as  second-order 
dynamical  systems: 

Cjs )  _  <*l  fB  ,, 

*(«>  '  s2  + 2Co>n+o2n  '  * 

where  C(s)  is  the  output  and  R(s)  is  the  input.  f  and 
are  referred  to  as  the  damping  ratio  and  the  natural 
frequency ,  respect i ve ly . 

The  values  for  f  and  ua  were  determined  by  experiment. 
The  control  vanes  were  mounted  on  the  servos  and  position 
data  was  obtained  for  a  range  of  step  input  commands.  The 
input  commands  ranged  from  three  to  thirty  degrees  of  vane 
deflection.  By  plotting  the  servo  angle  response  data 
versus  time,  a  series  of  step  input  response  curves  was 
developed.  Using  MATLAB,  trial  and  error  was  used  to  fit  a 
second  order  prototype  system  to  the  step  response  curves  of 
the  servos.  Figure  B.l  shows  the  "best  fit”  model.  This 
model  is  defined  by  f=0.707  and  w0=52.4  radians/second. 
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Thus,  the  second  order  servo  model  used  in  the  design  of  the 
Archytas  control  systems  is: 


C(S)  =  2745  ■  8 

R(s)  s2+74 . ls+2745 .8 


i 
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